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Abstract 

This paper is the outgrowth of lectures the author gave at the 
Physics Institute and the Chemistry Institute of the University of Sao 
Paulo at Sao Carlos, Brazil, and at the VHI'th Summer School on 
Electronic Structure of the Brazilian Physical Society. It is an attempt 
to introduce density-functional theory (DFT) in a language accessi- 
ble for students entering the field or researchers from other fields. It 
is not meant to be a scholarly review of DFT, but rather an infor- 
mal guide to its conceptual basis and some recent developments and 
advances. The Hohenberg-Kohn theorem and the Kohn-Sham equa- 
tions are discussed in some detail. Approximate density functionals, 
selected aspects of applications of DFT, and a variety of extensions 
of standard DFT are also discussed, albeit in less detail. Through- 
out it is attempted to provide a balanced treatment of aspects that 
are relevant for chemistry and aspects relevant for physics, but with 
a strong bias towards conceptual foundations. The paper is intended 
to be read before (or in parallel with) one of the many excellent more 
technical reviews available in the literature. 
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1 Preface 



This paper is the outgrowth of lectures the author gave at the Physics 
Institute and the Chemistry Institute of the University of Sao Paulo at Sao 
Carlos, Brazil, and at the VHI'th Summer School on Electronic Structure of 
the Brazilian Physical Society [I]. The main text is a description of density- 
functional theory (DFT) at a level that should be accessible for students 
entering the field or researchers from other fields. A large number of footnotes 
provides additional comments and explanations, often at a slightly higher 
level than the main text. A reader not familiar with DFT is advised to skip 
most of the footnotes, but a reader familiar with it may find some of them 
useful. 

The paper is not meant to be a scholarly review of DFT, but rather an 
informal guide to its conceptual basis and some recent developments and 
advances. The Hohenberg-Kohn theorem and the Kohn-Sham equations are 
discussed in some detail. Approximate density functionals, selected aspects 
of applications of DFT, and a variety of extensions of standard DFT are 
also discussed, albeit in less detail. Throughout it is attempted to provide 
a balanced treatment of aspects that are relevant for chemistry and aspects 
relevant for physics, but with a strong bias towards conceptual foundations. 
The text is intended to be read before (or in parallel with) one of the many 
excellent more technical reviews available in the literature. The author apol- 
ogizes to all researchers whose work has not received proper consideration. 
The limits of the author's knowledge, as well as the limits of the available 
space and the nature of the intended audience, have from the outset prohib- 
ited any attempt at comprehensiveness 

1 A first version of this text was published in 2002 as a chapter in the proceedings of 
the VlU'th Summer School on Electronic Structure of the Brazilian Physical Society [1]. 
The text was unexpectedly well received, and repeated requests from users prompted the 
author to electronically publish revised, updated and extended versions in the preprint 
archive \http:/ / arxiv. org/ archiv e/ cond-ma^ where the second (2003), third (2004) and 
fourth (2005) versions were deposited under the reference number |cond-mat /02 1 1443 
The present fifth (2006) version of this text, published in the Brazilian Journal of Physics, 
is approximately 50% longer than the first. Although during the consecutive revisions 
many embarrassing mistakes have been removed, and unclear passages improved upon, 
many other doubtlessly remain, and much beautiful and important work has not been 
mentioned even in passing. The return from electronic publishing to printed publishing, 
however, marks the completion of a cycle, and is intended to also mark the end of the 
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2 What is density-functional theory? 



Density-functional theory is one of the most popular and successful quan- 
tum mechanical approaches to matter. It is nowadays routinely applied for 
calculating, e.g., the binding energy of molecules in chemistry and the band 
structure of solids in physics. First applications relevant for fields tradition- 
ally considered more distant from quantum mechanics, such as biology and 
mineralogy are beginning to appear. Superconductivity, atoms in the focus 
of strong laser pulses, relativistic effects in heavy elements and in atomic nu- 
clei, classical liquids, and magnetic properties of alloys have all been studied 
with DFT. 

DFT owes this versatility to the generality of its fundamental concepts 
and the flexibility one has in implementing them. In spite of this flexibility 
and generality, DFT is based on quite a rigid conceptual framework. This 
section introduces some aspects of this framework in general terms. The 
following two sections, [3] and HI then deal in detail with two core elements 
of DFT, the Hohenberg-Kohn theorem and the Kohn-Sham equations. The 
final two sections, [5] and [6], contain a (necessarily less detailed) description of 
approximations typically made in practical DFT calculations, and of some 
extensions and generalizations of DFT. 

To get a first idea of what density-functional theory is about, it is useful to 
take a step back and recall some elementary quantum mechanics. In quantum 
mechanics we learn that all information we can possibly have about a given 
system is contained in the system's wave function, Here we will exclusively 
be concerned with the electronic structure of atoms, molecules and solids. 
The nuclear degrees of freedom (e.g., the crystal lattice in a solid) appear 
only in the form of a potential v(r) acting on the electrons, so that the wave 
function depends only on the electronic coordinates]! Nonrelativistically, this 
wave function is calculated from Schrodinger's equation, which for a single 
electron moving in a potential v(r) reads 



If there is more than one electron (i.e., one has a many-body problem) 

author's work on the Bird's-Eye View of Density- Functional Theory. 

2 This is the so-called Born-Oppenheimer approximation. It is common to call v(r) a 
'potential' although it is, strictly speaking, a potential energy. 




(1) 
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Schrodinger's equation becomes 



N ( h 2 v 2 \ 



*(n, r 2 . . . , r N ) = E^(r h r 2 . . . , r 



(2) 

where iV is the number of electrons and U(Ti,Tj) is the electron-electron 
interaction. For a Coulomb system (the only type of system we consider 
here) one has 

£ = = £iT37T- ( 3 ) 

i<j i<j \ Ll L n 

Note that this is the same operator for any system of particles interacting 
via the Coulomb interaction, just as the kinetic energy operator 

is the same for any nonrelativistic systemJl Whether our system is an atom, 
a molecule, or a solid thus depends only on the potential vfa). For an atom, 
e.g., 

fr = £*k) = £|r%T> (5) 

where Q is the nuclear charg^l and R the nuclear position. When dealing 
with a single atom, R is usually taken to be the zero of the coordinate system. 
For a molecule or a solid one has 

* = £«(**) = £]r%T. ( 6 ) 

where the sum on k extends over all nuclei in the system, each with charge 
Qk = Z^e and position R^. It is only the spatial arrangement of the R^ 
(together with the corresponding boundary conditions) that distinguishes, 
fundamentally, a molecule from a solid! Similarly, it is only through the term 



3 For materials containing atoms with large atomic number Z, accelerating the electrons 
to relativistic velocities, one must include relativistic effects by solving Dirac's equation 
or an approximation to it. In this case the kinetic energy operator takes a different form. 

In terms of the elementary charge e > and the atomic number Z, the nuclear charge 
is Q = Ze and the charge on the electron is q — —e. 

5 One sometimes says that T and U are 'universal', while V is system-dependent, or 
'nonuniversal'. We will come back to this terminology. 
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U that the (essentially simple) single-body quantum mechanics of Eq. (pQ) 
differs from the extremely complex many-body problem posed by Eq. (j2J). 
These properties are built into DFT in a very fundamental way. 

The usual quantum-mechanical approach to Schrodinger's equation (SE) 
can be summarized by the following sequence 



v(r) =^1> ^(ri, r 2 . . . , r N ) ^^=k4>^ observables, 



(7) 



i.e., one specifies the system by choosing v(r), plugs it into Schrodinger's 
equation, solves that equation for the wave function and then calculates 
observables by taking expectation values of operators with this wave function. 
One among the observables that are calculated in this way is the particle 
density 



n(r) = N 



J d 3 r 2 J d 3 r 3 ...J d 3 r N V*(r, r 2 . . . , r N )V(r, r 2 ...,r N ). (8) 



Many powerful methods for solving Schrodinger's equation have been devel- 
oped during decades of struggling with the many-body problem. In physics, 
for example, one has diagrammatic perturbation theory (based on Feynman 
diagrams and Green's functions), while in chemistry one often uses config- 
uration interaction (CI) methods, which are based on systematic expansion 
in Slater determinants. A host of more special techniques also exists. The 
problem with these methods is the great demand they place on one's compu- 
tational resources: it is simply impossible to apply them efficiently to large 
and complex systems. Nobody has ever calculated the chemical properties 
of a 100-atom molecule with full CI, or the electronic structure of a real 
semiconductor using nothing but Green's functions]^ 



6 A simple estimate of the computational complexity of this task is to imagine a real- 
space representation of $ on a mesh, in which each coordinate is discretized by using 
20 mesh points (which is not very much). For N electrons, 4" becomes a function of 3N 
coordinates (ignoring spin, and taking *5> to be real), and 20 3N values are required to 
describe ^ on the mesh. The density n(r) is a function of three coordinates, and requires 
20 3 values on the same mesh. CI and the Kohn-Sham formulation of DFT additionally 
employ sets of single-particle orbitals. N such orbitals, used to build the density require 
20 3 iV values on the same mesh. (A CI calculation employs also unoccupied orbitals, and 
requires more values.) For N = 10 electrons, the many-body wave function thus requires 
20 30 /20 3 « 10 35 times more storage space than the density, and 20 30 /(10 x 20 3 ) « 10 34 
times more than sets of single-particle orbitals. Clever use of symmetries can reduce these 
ratios, but the full many-body wave function remains unaccessible for real systems with 
more than a few electrons. 
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It is here where DFT provides a viable alternative, less accurate per- 
hapsJ3 but much more versatile. DFT explicitly recognizes that nonrela- 
tivistic Coulomb systems differ only by their potential v(r), and supplies a 
prescription for dealing with the universal operators T and U once and for 
allH Furthermore, DFT provides a way to systematically map the many- 
body problem, with U, onto a single-body problem, without U. All this is 
done by promoting the particle density n(r) from just one among many ob- 
servables to the status of key variable, on which the calculation of all other 
observables can be based. This approach forms the basis of the large ma- 
jority of electronic-structure calculations in physics and chemistry. Much of 
what we know about the electrical, magnetic, and structural properties of 
materials has been calculated using DFT, and the extent to which DFT has 
contributed to the science of molecules is reflected by the 1998 Nobel Prize 
in Chemistry, which was awarded to Walter Kohn [3], the founding father 
of DFT, and John Pople [I], who was instrumental in implementing DFT in 
computational chemistry. 

The density-functional approach can be summarized by the sequence 



n(r) =^ ^(n,...,r N ) =>>v(r), (9) 



i.e., knowledge of n(r) implies knowledge of the wave function and the po- 
tential, and hence of all other observables. Although this sequence describes 
the conceptual structure of DFT, it does not really represent what is done in 
actual applications of it, which typically proceed along rather different lines, 
and do not make explicit use of many-body wave functions. The following 
chapters attempt to explain both the conceptual structure and some of the 



7 Accuracy is a relative term. As a theory, DFT is formally exact. Its performance 
in actual applications depends on the quality of the approximate density functionals em- 
ployed. For small numbers of particles, or systems with special symmetries, essentially 
exact solutions of Schrodinger's equation can be obtained, and no approximate functional 
can compete with exact solutions. For more realistic systems, modern (2005) sophisticated 
density functionals attain rather high accuracy. Data on atoms are collected in Table [1] 
in Sec. 15.21 Bond-lengths of molecules can be predicted with an average error of less than 
O.OOlnm, lattice constants of solids with an average error of less than 0.005nm, and molec- 
ular energies to within less than 0.2eV [2]. (For comparison: already a small molecule, 
such as water, has a total energy of 2081. leV). On the other hand, energy gaps in solids 
can be wrong by 100%! 

8 We will see that in practice this prescription can be implemented only approximately. 
Still, these approximations retain a high degree of universality in the sense that they often 
work well for more than one type of system. 
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many possible shapes and disguises under which this structure appears in 
applications. 

The literature on DFT is large, and rich in excellent reviews and overviews. 
Some representative examples of full reviews and systematic collections of 
research papers are Refs. [5-19]. The present overview of DFT is much less 
detailed and advanced than these treatments. Introductions to DFT that 
are more similar in spirit to the present one (but differ in emphasis and se- 
lection of topics) are the contribution of Levy in Ref. [9], the one of Kurth 
and Per dew in Refs. [15] and [16] , and Ref. [20] by Makov and Argaman. My 
aim in the present text is to give a bird's-eye view of DFT in a language that 
should be accessible to an advanced undergraduate student who has com- 
pleted a first course in quantum mechanics, in either chemistry or physics. 
Many interesting details, proofs of theorems, illustrative applications, and 
exciting developments had to be left out, just as any discussion of issues that 
are specific to only certain subfields of either physics or chemistry. All of 
this, and much more, can be found in the references cited above, to which 
the present little text may perhaps serve as a prelude. 

3 DFT as a many-body theory 



3.1 Functionals and their derivatives 

Before we discuss density-functional theory more carefully, let us introduce 
a useful mathematical tool. Since according to the above sequence the wave 
function is determined by the density, we can write it as ^ = ^[n]^, r 2 , . . . r N ), 
which indicates that ^ is a function of its iV spatial variables, but a functional 
of n(r). 

Functionals. More generally, a functional F[n] can be defined (in an 
admittedly mathematically sloppy way) as a rule for going from a function 
to a number, just as a function y = f(x) is a rule (/) for going from a number 
(x) to a number (y). A simple example of a functional is the particle number, 



which is a rule for obtaining the number N, given the function n(r). Note 
that the name given to the argument of n is completely irrelevant, since the 




(10) 
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functional depends on the function itself, not on its variable. Hence we do 
not need to distinguish _F[n(r)] from, e.g., F[n(r')]. Another important case 
is that in which the functional depends on a parameter, such as in 




(11) 



which is a rule that for any value of the parameter r associates a value 
■Utf[n](r) with the function n(r'). This term is the so-called Hartree potential, 
which we will repeatedly encounter below. 

Functional variation. Given a function of one variable, y = f{x), one can 
think of two types of variations of y, one associated with x, the other with /. 
For a fixed functional dependence f(x), the ordinary differential dy measures 
how y chang result of a variation x — > x + dx of the variable x. This 

is the variation studied in ordinary calculus. Similarly, for a fixed point x, 
the functional variation Sy measures how the value y at this point changes 
as a result of a variation in the functional form f(x). This is the variation 
studied in variational calculus. 

Functional derivative. The derivative formed in terms of the ordinary 
differential, df /dx, measures the first-order change of y — f(x) upon changes 
of x, i.e., the slope of the function f(x) at x: 



The functional derivative measures, similarly, the first-order change in a func- 
tional upon a functional variation of its argument: 



where the integral arises because the variation in the functional F is deter- 
mined by variations in the function at all points in space. The first-order 
coefficient (or 'functional slope') s(x) is defined to be the functional derivative 



The functional derivative allows us to study how a functional changes 
upon changes in the form of the function it depends on. Detailed rules for 
calculating functional derivatives are described in Appendix A of Ref. [6]. A 
general expression for obtaining functional derivatives with respect to n(x) of 
a functional F[n] = j f(n, n', n", n'", x)dx, where primes indicate ordinary 




(12) 



F[f(x) + 6f{x)\ = F[f{x)\ + J s{x) 5f(x) dx + 0(5 f), (13) 



5F[f]/5f(x). 
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derivatives of n{x) with respect to x, is [5] 

SF[n] = d£_ + d?_ df_ _ d^d£_ + 

Sn(x) dn dx dn' dx 2 dn" dx 3 dn'" 

This expression is frequently used in DFT to obtain xc potentials from xc 
energies @ 



3.2 The Hohenberg-Kohn theorem 

At the heart of DFT is the Hohenberg-Kohn (HK) theorem. This theo- 
rem states that for ground states Eq. (jSJ) can be inverted: given a ground- 
state density n (r) it is possible, in principle, to calculate the corresponding 
ground-state wave function *o( r i) r 2 • • • , y n)- This means that ^/q is a func- 
tional of n . Consequently, all ground-state observables are functionals of 
no, too. If \l/o can be calculated from n and vice versa, both functions are 
equivalent and contain exactly the same information. At first sight this seems 
impossible: how can a function of one (vectorial) variable r be equivalent to a 
function of N (vectorial) variables ri . . .rjv? How can one arbitrary variable 
contain the same information as N arbitrary variables? 

The crucial fact which makes this possible is that knowledge of no(r) 
implies implicit knowledge of much more than that of an arbitrary func- 
tion f(r). The ground-state wave function \l/ must not only reproduce the 
ground-state density, but also minimize the energy. For a given ground-state 
density n (r), we can write this requirement as 

E v0 = min(*|f + U + VI*), (15) 

where E Vt0 denotes the ground-state energy in potential v(r). The preceding 
equation tells us that for a given density n (r) the ground-state wave function 
\l/o is that which reproduces this rio(r) and minimizes the energy. 
For an arbitrary density n(r), we define the functional 

E v [n] = min(*|f + U + V\ty). (16) 

9 The use of functionals and their derivatives is not limited to density-functional 
theory, or even to quantum mechanics. In classical mechanics, e.g., one expresses the 
Lagrangian C in terms of of generalized coordinates q(x, t) and their temporal deriva- 
tives q(x,t), and obtains the equations of motion from extremizing the action functional 
A[q] — J C(q, q; t)dt. The resulting equations of motion are the well-known Euler-Lagrange 
equations = j^M = |p — ^ ^ , which are a special case of Eq. (TTJ)) . 
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If n is a density different from the ground-state density n in potential v(r), 
then the \I/ that produce this n are different from the ground-state wave 
function ^o, and according to the variational principle the minimum obtained 
from E v [n] is higher than (or equal to) the ground-state energy E Vt0 = E v [n }. 
Thus, the functional E v [n] is minimized by the ground-state density n , and 
its value at the minimum is E v $. 

The total-energy functional can be written as 



E v [n] = min^ n (*|f + U\V) + J d 3 r n(r)v(r) =: F[n] + V[n], (17) 



where the internal-energy functional F[n] = min^^ n (\l/|T + U\^) is inde- 
pendent of the potential f(r), and thus determined only by the structure of 
the operators U and T. This universality of the internal-energy functional 
allows us to define the ground-state wave function \l/ as that antisymmetric 
iV-particle function that delivers the minimum of F[n] and reproduces no- 
If the ground state is nondegenerate (for the case of degeneracy see foot- 
note 12), this double requirement uniquely determines \l/o in terms of no(r), 
without having to specify w(r) explicitly] 10 ! 

Equations (|T5l) to (|T71) constitute the constrained-search proof of the 
Hohenberg-Kohn theorem, given independently by M. Levy [22] and E. Lieb 
|23j . The original proof by Hohenberg and Kohn [23] proceeded by assuming 
that \l/o was n °t determined uniquely by n and showed that this produced a 
contradiction to the variational principle. Both proofs, by constrained search 
and by contradiction, are elegant and simple. In fact, it is a bit surprising 
that it took 38 years from Schrodinger's first papers on quantum mechanics 
[25] to Hohenberg and Kohn's 1964 paper containing their famous theorem 



Since 1964, the HK theorem has been thoroughly scrutinized, and several 
alternative proofs have been found. One of these is the so-called 'strong form 
of the Hohenberg-Kohn theorem', based on the inequality [26 ] 1271 [28] 

d 3 rAn(r)Av (r) < 0. (18) 

Here Av(r) is a change in the potential, and An(r) is the resulting change in 
the density. We see immediately that if Av ^ we cannot have An(r) = 0, 



10 Note that this is exactly the opposite of the conventional prescription to specify the 
Hamiltonian via v(r), and obtain from solving Schrodinger's equation, without having 
to specify n(r) explicitly. 
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i.e., a change in the potential must also change the density. This observation 
implies again the HK theorem for a single density variable: there cannot 
be two local potentials with the same ground-state charge density. A given 
iV-particle ground-state density thus determines uniquely the corresponding 
potential, and hence also the wave function. Moreover, (|18p establishes a 
relation between the signs of An(r) and Av(r): if Av is mostly positive, 
An(r) must be mostly negative, so that their integral over all space is nega- 
tive. This additional information is not immediately available from the two 
classic proofs, and is the reason why this is called the 'strong' form of the 
HK theorem. Equation (JT81) can be obtained along the lines of the standard 
HK proof [251 EZ] , but it can be turned into an independent proof of the HK 
theorem because it can also be derived perturbatively (see, e.g., section 10.10 
of Ref. [28]). 

Another alternative argument is valid only for Coulomb potentials. It is 
based on Kato's theorem, which states [29], [30] that for such potentials the 
electron density has a cusp at the position of the nuclei, where it satisfies 



ao dn(r) 



2n(r) dr 



(19) 

+R fc 



Here R& denotes the positions of the nuclei, their atomic number, and 
a o — h 2 /me 2 is the Bohr radius. For a Coulomb system one can thus, in 
principle, read off all information necessary for completely specifying the 
Hamiltonian directly from examining the density distribution: the integral 
over n (r) yields N, the total particle number; the position of the cusps of 
n(r) are the positions of the nuclei, R&; and the derivative of n(r) at these 
positions yields Zk by means of Eq. ffT9l) . This is all one needs to specify the 
complete Hamiltonian of Eq. ([2]) (and thus implicitly all its eigenstates). In 
practice one almost never knows the density distribution sufficiently well to 
implement the search for the cusps and calculate the local derivatives. Still, 
Kato's theorem provides a vivid illustration of how the density can indeed 
contain sufficient information to completely specify a nontrivial Hamilto- 
nian 

For future reference we now provide a commented summary of the content 
of the HK theorem. This summary consists of four statements: 



11 Note that, unlike the full Hohcnbcrg-Kohn theorem, Kato's theorem does apply only 
to superpositions of Coulomb potentials, and can therefore not be applied directly to the 
effective Kohn-Sham potential. 
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(1) The nondegenerate ground-state (GS) wave function is a unique func- 
tional of the GS density o 



^ {r u r 2 ...,r N ) = #[n (r)]. 



(20) 



This is the essence of the HK theorem. As a consequence, the GS expectation 
value of any observable O is a functional of no(r), too: 



O = 0[no) = (*[7io]|0|*[no]>. 



(21) 



(2) Perhaps the most important observable is the GS energy. This energy 
E vfi = E v [n ] = (tt[no]|#|*[no]), (22) 
where H = T + U + V, has the variational propertjQ 

(23) 



E v [n ) < E v [n'} 



where n is GS density in potential V and n' is some other density. This 
is very similar to the usual variational principle for wave functions. From 
a calculation of the expectation value of a Hamiltonian with a trial wave 
function \P' that is not its GS wave function ty one can never obtain an 
energy below the true GS energy, 



E Vj0 = E V [V ] = (*o|£|*o> < (V'\H\V') = E v [¥\. 



(24) 



Similarly, in exact DFT, if E[n] for fixed v ext is evaluated for a density that 
is not the GS density of the system in potential v ext , one never finds a result 
below the true GS energy. This is what Eq. fl23|) says, and it is so impor- 
tant for practical applications of DFT that it is sometimes called the second 
Hohenberg-Kohn theorem (Eq. f[2~Tl) is the first one, then). 



12 If the ground state is degenerate, several of the degenerate ground-state wave functions 
may produce the same density, so that a unique functional ^S[n] does not exist, but by 
definition these wave functions all yield the same energy, so that the functional i£u[n] 
continues to exist and to be minimized by no. A universal functional F[n] can also still 
be defined [5]. 

13 The minimum of E[n] is thus attained for the ground-state density. All other extrema 
of this functional correspond to densities of excited states, but the excited states obtained 
in this way do not necessarily cover the entire spectrum of the many-body Hamiltonian 
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In performing the minimization of E v [n] the constraint that the total par- 
ticle number N is an integer is taken into account by means of a Lagrange 
multiplier, replacing the constrained minimization of E v [n] by an uncon- 
strained one of E v [n] — iiN. Since N — J d 3 rn(r), this leads to 



SE v [n] 
5n(r) 



dE 



(25) 



where 11 is the chemical potential. 

(3) Recalling that the kinetic and interaction energies of a nonrelativistic 
Coulomb system are described by universal operators, we can also write E v 

as 

(26) 



E v [n 


= T 


n 


+ U[n] + V[n] = F[n] + V[n], 



where T[n] and U[n] are universal Junctionals [defined as expectation values 
of the type (l2Tj) of T and U], independent of v(r). On the other hand, the 
potential energy in a given potential v(r) is the expectation value of Eq. (jBJ), 



V[n] = / a?r n(r)v(r), 



(27) 



and obviously nonuniversal (it depends on v(r), i.e., on the system under 
study), but very simple: once the system is specified, i.e., v(r) is known, the 
functional V[n] is known explicitly. 

(4) There is a fourth substatement to the HK theorem, which shows 
that if v(r) is not hold fixed, the functional V[n] becomes universal: the GS 
density determines not only the GS wave function \l/ 0) but, up to an additive 
constant, also the potential V = V[n }. This is simply proven by writing 
Schrodinger's equation as 



V(Ti 



Eu — 



(28) 



which shows that any eigenstate (and thus in particular the ground state 
= ^[^o]) determines the potential operator V up to an additive constant, 
the corresponding eigenenergy. As a consequence, the explicit reference to 
the potential v in the energy functional E v [n] is not necessary, and one can 
rewrite the ground-state energy as 



E = E[n ] = mn Q ]\T + U + V[n }mn }). 



(29) 
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Another consequence is that no now does determine not only the GS wave 
function but the complete Hamiltonian (the operators T and U are fixed), 
and thus all excited states, too: 



* fe (ri,r 2 .. -,r N ) = * fc [n ], 



(30) 



where k labels the entire spectrum of the many-body Hamiltonian H . 



3.3 Complications: N and ^-representability of densi- 
ties, and nonuniqueness of potentials 

Originally the fourth statement was considered to be as sound as the other 
three. However, it has become clear very recently, as a consequence of work of 
H. Eschrig and W. Pickett |32j and, independently, of the author with G. Vi- 
gnale [331 IS], that there are significant exceptions to it. In fact, the fourth 
substatement holds only when one formulates DFT exclusively in terms of the 
charge density, as we have done up to this point. It does not hold when one 
works with spin densities (spin-DFT) or current densities (current-DFT) f^l 
In these (and some other) cases the densities still determine the wave func- 
tion, but they do not uniquely determine the corresponding potentials. This 
so-called nonuniqueness problem has been discovered only recently, and its 
consequences are now beginning to be explored J27J E2l [331 IM1 [351 USE [37J [38] ■ 
It is clear, however, that the fourth substatement is, from a practical point 
of view, the least important of the four, and most applications of DFT do 
not have to be reconsidered as a consequence of its eventual failure. (But 
some do: see Refs. [331 [M] f° r examples.) 

Another conceptual problem with the HK theorem, much better known 
and more studied than nonuniqueness, is representability. To understand 
what representability is about, consider the following two questions: (i) How 
does one know, given an arbitrary function n(r), that this function can be 
represented in the form (jSj), i.e., that it is a density arising from an antisym- 
metric A-body wave function \l/(ri . . . rjv)? (ii) How does one know, given a 
function that can be written in the form (jSJ), that this density is a ground- 
state density of a local potential u(r)? The first of these questions is known as 
the A-representability problem, the second is called v-representability. Note 
that these are quite important questions: if one should find, for example, in a 

14 In Section [5] we will briefly discuss these formulations of DFT. 
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numerical calculation, a minimum of E v [n] that is not iV-representable, then 
this minimum is not the physically acceptable solution to the problem at 
hand. Luckily, the iV-representability problem of the single-particle density 
has been solved: any nonnegative function can be written in terms of some 
antisymmetric \&(ri, r 2 . . . , r/v) in the form (|8]) [391 ED]- 

No similarly general solution is known for the w-representability problem. 
(The HK theorem only guarantees that there cannot be more than one po- 
tential for each density, but does not exclude the possibility that there is less 
than one, i.e., zero, potentials capable of producing that density.) It is known 
that in discretized systems every density is ensemble t>-representable, which 
means that a local potential with a degenerate ground state can always be 
found, such that the density n(r) can be written as linear combination of 
the densities arising from each of the degenerate ground states [HI H2J H3] . 
It is not known if one of the two restrictions ('discretized systems', and 
'ensemble') can be relaxed in general (yielding 'in continuum systems' and 
'pure-state' respectively), but it is known that one may not relax both: there 
are densities in continuum systems that are not representable by a single non- 
degenerate antisymmetric function that is ground state of a local potential 
v(t) [SI SH E2J H3] . In any case, the constrained search algorithm of Levy and 
Lieb shows that f-representability in an interacting system is not required 
for the proof of the HK theorem. For the related question of simultaneous 
f-representability in a noninteracting system, which appears in the context 
of the Kohn-Sham formulation of DFT, see footnotes 34 and 35. 

3.4 A preview of practical DFT 

After these abstract considerations let us now consider one way in which one 
can make practical use of DFT. Assume we have specified our system (i.e., 
v(r) is known). Assume further that we have reliable approximations for 
U[n] and T[n]. In principle, all one has to do then is to minimize the sum of 
kinetic, interaction and potential energies 

E v [n] = T[n] + U[n) + V[n] = T[n) + U[n) + J d 3 rn(r)v(r) (31) 

with respect to n(r). The minimizing function no(r) is the system's GS 
charge density and the value E Vt0 = E v [n ] is the GS energy. Assume now 
that v(r) depends on a parameter a. This can be, for example, the lattice 
constant in a solid or the angle between two atoms in a molecule. Calculation 
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of E v<0 for many values of a allows one to plot the curve E Vt o(a) and to find 
the value of a that minimizes it. This value, do, is the GS lattice constant or 
angle. In this way one can calculate quantities like molecular geometries and 
sizes, lattice constants, unit cell volumes, charge distributions, total energies, 
etc. By looking at the change of £^0(0) with a one can, moreover, calculate 
compressibilities, phonon spectra and bulk moduli (in solids) and vibrational 
frequencies (in molecules). By comparing the total energy of a composite 
system (e.g., a molecule) with that of its constituent systems (e.g., individual 
atoms) one obtains dissociation energies. By calculating the total energy 
for systems with one electron more or less one obtains electron affinities 
and ionization energies!^! By appealing to the Hellman-Feynman theorem 
one can calculate forces on atoms from the derivative of the total energy 
with respect to the nuclear coordinates. All this follows from DFT without 
having to solve the many-body Schrodinger equation and without having 
to make a single-particle approximation. For brief comments on the most 
useful additional possibility to also calculate single-particle band structures 
see Sees, and 

In theory it should be possible to calculate all observables, since the 
HK theorem guarantees that they are all functionals of no(r). In practice, 
one does not know how to do this explicitly. Another problem is that the 
minimization of E v [n] is, in general, a tough numerical problem on its own. 
And, moreover, one needs reliable approximations for T[n] and U[n] to begin 
with. In the next section, on the Kohn-Sham equations, we will see one widely 
used method for solving these problems. Before looking at that, however, it 
is worthwhile to recall an older, but still occasionally useful, alternative: the 
Thomas- Fermi approximation. 

In this approximation one sets 

U[n] « U H [n] = ^jd 3 rj aV ^^ , (32) 

i.e., approximates the full interaction energy by the Hartree energy, the elec- 
trostatic interaction energy of the charge distribution n(r). One further 

Electron affinities are typically harder to obtain than ionization energies, because 
within the local-density and generalized-gradient approximations the N + l'st electron is 
too weakly bound or even unbound: the asymptotic effective potential obtained from these 
approximations decays exponentially, and not as 1/r, i.e., it approaches zero so fast that 
binding of negative ions is strongly suppressed. Self-interaction corrections or other fully 
nonlocal functionals are needed to improve this behaviour. 
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approximates, initially, 

T[n] « T LDA [n] = J d 3 rt hom (n(r)), (33) 

where t hom (n) is the kinetic-energy density of a homogeneous interacting 
system with (constant) density n. Since it refers to interacting electrons 
t hom (n) is not known explicitly, and further approximations are called for. 
As it stands, however, this formula is already a first example of a local- 
density approximation (LDA). In this type of approximation one imagines 
the real inhomogeneous system (with density n(r) in potential v(r)) to be 
decomposed in small cells in each of which n(r) and v(r) are approximately 
constant. In each cell (i.e., locally) one can then use the per-volume energy 
of a homogeneous system to approximate the contribution of the cell to the 
real inhomogeneous one. Making the cells infinitesimally small and summing 
over all of them yields Eq. (I33p . 

For a noninteracting system (specified by subscript s, for 'single-particle') 
the function t h s om {n) is known explicitly, t h s om {n) = 3ft 2 (3vr 2 ) 2 / 3 n 5 / 3 /(10m) 
(see also Sec. 15.11) . This is exploited to further approximate 

T[n\ « T LDA [n] « T s LDA [n] = J d 3 rt h s om (n(r)), (34) 

where T S LDA [ n] is the local-density approximation to T a [n], the kinetic energy 
of noninteracting electrons of density n. Equivalently, it may be considered 
the noninteracting version of T LDA [n}. (The quantity T s [n] will reappear 
below, in discussing the Kohn-Sham equations.) The Thomas- Fermi approx- 
imation^ then consists in combining (1321) and (154]) and writing 

E[n] = T[n] + U[n] + V[n] « E TF [n] = T s LDA [n] + U H [n] + V[n}. (35) 

A major defect of the Thomas- Fermi approximation is that within it molecules 
are unstable: the energy of a set of isolated atoms is lower than that of the 
bound molecule. This fundamental deficiency, and the lack of accuracy result- 
ing from neglect of correlations in (|32|) and from using the local approxima- 
tion (1341) for the kinetic energy, limit the practical use of the Thomas-Fermi 

16 The Thomas-Fermi approximation for screening, discussed in many books on solid- 
state physics, is obtained by minimizing E TF [n] with respect to n and linearizing the 
resulting relation between v(r) and n(r). It thus involves one more approximation (the 
linearization) compared to what is called the Thomas-Fermi approximation in DFT |44|. 
In two dimensions no linearization is required and both become equivalent [44| . 
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approximation in its own right. However, it is found a most useful starting 
point for a large body of work on improved approximations in chemistry and 
physics [12, 30J. More recent approximations for T[n] can be found, e.g., in 
Refs. |l5l H6l H7], in the context of orbital- free DFT. The extension of the 
local-density concept to the exchange-correlation energy is at the heart of 
many modern density functionals (see Sec. 15.11) . 

3.5 From wave functions to density functionals via Green's 
functions and density matrices 

It is a fundamental postulate of quantum mechanics that the wave function 
contains all possible information about a system in a pure state at zero tem- 
perature, whereas at nonzero temperature this information is contained in 
the density matrix of quantum statistical mechanics. Normally, this is much 
more information that one can handle: for a system with iV = 100 parti- 
cles the many-body wave function is an extremely complicated function of 
300 spatial and 100 spirit variables that would be impossible to manipulate 
algebraically or to extract any information from, even if it were possible to 
calculate it in the first place. For this reason one searches for less complicated 
objects to formulate the theory. Such objects should contain the experimen- 
tally relevant information, such as energies, densities, etc., but do not need 
to contain explicit information about the coordinates of every single parti- 
cle. One class of such objects are Green's functions, which are described in 
the next subsection, and another are reduced density matrices, described in 
the subsection 13.5.21 Their relation to the wave function and the density is 
summarized in Fig. [IJ 

3.5.1 Green's functions 

Readers with no prior knowledge of (or no interest in) Green's functions 
should skip this subsection, which is not necessary for understanding the 
following sections. 

In mathematics one usually defines the Green's function of a linear op- 
erator C via [z — C(r))G(x,x';z) = S(x — x'), where 5(x — x') is Dirac's 

17 To keep the notation simple, spin labels are either ignored or condensed into a common 
variable x :— (rs) in most of this text. They will only be put back explicitly in discussing 
spin-density-functional theory, in Sec. [BJ 
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delta function. For a single quantum particle in potential v(r) one has, for 
example, 

T, 2 X7 2 1 

G(°)(r,r';£) = M(r-r'). (36) 



E 



2m 



v(r) 



Many applications of such single-particle Green's functions are discussed in 
Ref. [21]. In many-body physics it is useful to also introduce more compli- 
cated Green's functions. In an interacting system the single-particle Green's 
function is modified by the presence of the interaction between the particles 1^1 
In general it now satisfies the equation^! 



dt + ~2m~ ~ V{r) 



G{r,t; r',t') = H5(r - r')5(t - t') 
J d 3 xU(r -x)G {2) (rt,xt; r't',xt+). 



(37) 
Only 



where G^ 2 \rt, xt; r't', xt + ) is the two-particle Green's function 
for noninteracting systems (U = 0) is G(r, t; r', t') a Green's function in the 
mathematical sense of the word. In terms of G(r,t;r',t') one can explicitly 
express the expectation value of any single-body operator (such as the po- 
tential, the kinetic energy, the particle density, etc.), and also that of certain 
two-particle operators, such as the Hamiltonian in the presence of particle- 
particle interactions] 20 ! 

One way to obtain the single-particle Green's function is via solution of 
what is called Dyson's equation [2T | B8l H9] . 

G(r,t;r r ,t') = G (0) (r,t; r', t') 
J d 3 x J d 3 x' J d 3 r J d 3 r' G (0) (r, £; x, r)S(x, r, x', r / )G(x / , r'; r', t'), (38) 



+ 



where S is known as the irreducible self energy [HI SSI SI] and G (0) is the 
Green's function in the absence of any interaction. This equation (which 



18 Note that expressions like 'two-particle operator' and 'single-particle Green's function' 
refer to the number of particles involved in the definition of the operator (two in the case 
of an interaction, one for a potential energy, etc.), not to the total number of particles 
present in the system. 

19 When energy is conserved, i.e., the Hamiltonian does not depend on time, G(r, t; r', t 1 ) 
depends on time only via the difference t—t' and can be written as G(r, r'; t—t'). By Fourier 
transformation with respect to t — t 1 one then passes from G(r, r'; t — t') to G(r, r'; E) of 
Eq. (J36]). 

20 For arbitrary two-particle operators one needs the full two-particle Green's function 
G( 2 ). 
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we will not attempt to solve here) has a characteristic property that we will 
meet again when we study the (much simpler) Kohn-Sham and Hartree-Fock 
equations, in Sec. HI the integral on the right-hand side, which determines 
G on the left-hand side, depends on G itself. The mathematical problem 
posed by this equation is thus nonlinear. We will return to such nonlinearity 
when we discuss self-consistent solution of the Kohn-Sham equation. The 
quantity £ will appear again in Sec. 14.2.31 when we discuss the meaning of 
the eigenvalues of the Kohn-Sham equation. 

The single-particle Green's function is related to the irreducible self en- 
ergy by Dyson's equation (138|) and to the two-particle Green's function by 
the equation of motion (1371) . It can also be related to the xc potential of 
DFT by the Sham-Schluter equation [SU] 

J d 3 r'v xc (r) J uG s (r,r'; u)G(r', r; to) = 
J d 3 r> J d 3 r"J uG s (r,r';u)E xc {r',r";u J )G{yr",r;u;), (39) 

where G s is the Green's function of noninteracting particles with density 
n(r) (i.e., the Green's function of the Kohn-Sham equation, see Sec. 14. 2j) . 
and Y* xc (r', r"; to) = S(r', r"; uj) — 5(r' — r")v H (r') represents all contributions 
to the full irreducible self energy beyond the Hartree potential. 

A proper discussion of £ and G requires a formalism known as second 
quantization [2TI HH] and usually proceeds via introduction of Feynman dia- 
grams. These developments are beyond the scope of the present overview. A 
related concept, density matrices, on the other hand, can be discussed easily. 
The next section is devoted to a brief description of some important density 
matrices. 



3.5.2 Density matrices 

For a general quantum system at temperature T, the density operator in a 
canonical ensemble is defined as 

7 = 6XP ^ , (40) 
Tr[exp-^] V ; 

where Tr[-} is the trace and (3 = l/{k B T). Standard textbooks on statistical 
physics show how this operator is obtained in other ensembles, and how it 
is used to calculate thermal and quantum expectation values. Here we focus 
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on the relation to density-functional theory. To this end we write 7 in the 
energy representation as 



where is eigenfunction of H, and the sum is over the entire spectrum of 
the system, each state being weighted by its Boltzmann weight exp~^ Ei . At 
zero temperature only the ground-state contributes to the sums, so that 

7=|**>W (42) 

The coordinate-space matrix element of this operator for an iV-particle sys- 
tem is 

(Xi,X 2 , ..X N \j\x[, x' 2 , ..x' N ) = ^{Xi, X 2 , ..X A r)*^(x / 1 , x' 2 , ..X N ) 

=: 7(2:1, x 2 , ..x N ; x[,x 2 , ..x' N ), (43) 

which shows the connection between the density matrix and the wave func- 
tion. (We use the usual abbreviation x = rs for space and spin coordi- 
nates.) The expectation value of a general iV-particle operator O is obtained 
from O = (O) = f dxi f dx 2 ... f dxN^(xi,x 2 , ..xn)*0^(xi,x 2 , ..xn), which 
for multiplicative operators becomes 

(O) = J dxi J dx 2 ... J dx N Oj(xi,x 2 ,..x N ;x 1 ,x 2 ,..x N ) (44) 

and involves only the function 7(21, x 2 , ..xn', x%, x 2 , ..xn), which is the di- 
agonal element of the matrix 7. Most operators we encounter in quantum 
mechanics are one or two-particle operators and can be calculated from re- 
duced density matrices, that depend on less than 2N variables^ The reduced 
two-particle density matrix is defined as 

72(xi,x 2 ;x / 1 ,x / 2 ) = 
' f dx 3 [ dx 4 ... [ dxN , y(xi, x 2 , £ 3 , x 4 , ..xn', x[, x' 2} x 3 , x 4 , ..xn), (45) 



21 Just as for Green's functions, expressions like 'two-particle operator' and 'two-particle 
density matrix' refer to the number of particles involved in the definition of the operator 
(two in the case of an interaction, one for a potential energy, etc.), not to the total number 
of particles present in the system. 
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where N(N — l)/2 is a convenient normalization factor. This density ma- 
trix determines the expectation value of the particle-particle interaction, 
of static correlation and response functions, of the xc hole, and some re- 
lated quantities. The pair-correlation function g(x,x'), e.g., is obtained from 
the diagonal element of 72(^1, x 2 ; x[, x' 2 ) according to 72(^1, x 2 ; x 2 ) =: 
n(xi)n(x 2 )g(x, x'). 

Similarly, the single-particle density matrix is defined as 

7(^1,^1) = 

iV J dx 2 J dx 3 J dx±... J dxN"f(xi, x 2 , x 3 , x 4 , ..xn', x[, x 2 , x 3 , x 4 , ..xn) 
= N J dx 2 J dx 3 J dx 4 ... J dx N ^*(x 1 ,x 2 ,x 3 , .., x N )^(x[, x 2 , x 3 , ..,x N ). (46) 

The structure of reduced density matrices is quite simple: all coordinates that 
7 does not depend upon are set equal in ^ and \f r *, and integrated over. The 
single-particle density matrix can also be considered the time-independent 
form of the single-particle Green's function, since it can alternatively be 
obtained from 

1 (x,x) = -ihuiG(x,x';t-t f ). (47) 

In the special case that the wave function ^ is a Slater determinant, i.e., 
the wave function of iV noninteracting fermions, the single-particle density 
matrix can be written in terms of the orbitals comprising the determinant, 

as 

7 (*y)=£0*(a;)<MA (48) 

3 

which is known as the Dirac (or Dirac-Fock) density matrix. 

The usefulness of the single-particle density matrix becomes apparent 
when we consider how one would calculate the expectation value of a mul- 
tiplicative single-particle operator A = X)f a(rj) (such as the potential V = 
EfVr,)): 



(A) = J dxi . . .J dx N ^*(xi,x 2 , -,xn) 



■ N 

^2a(xi) V(x 1 ,x 2 ,..,x N ) (49) 



= N J dxi . . . J dx N 1 &*(x 1 , x 2 , .., x N )a(x±)^(xi,x 2 , ..,x N ) (50) 

= J dxa(x) n /(x,x), (51) 
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which is a special case of Eq. (j44p . The second line follows from the first by 
exploiting that the fermionic wave function \l/ changes sign upon interchange 
of two of its arguments. The last equation implies that if one knows j(x,x) 
one can calculate the expectation value of any multiplicative single-particle 
operator in terms of it, regardless of the number of particles present in the 
systemic The simplification is enormous, and reduced density matrices are 
very popular in, e.g., computational chemistry for precisely this reason. More 
details are given in, e.g., Ref. [6]. The full density operator, Eq. (|40j) . on the 
other hand, is the central quantity of quantum statistical mechanics. 

It is not possible to express expectation values of two-particle operators, 
such as the interaction itself, or the full Hamiltonian (i.e., the total energy), 
explicitly in terms of the single-particle density matrix 7(1*, r'). For this 
purpose one requires the two-particle density matrix. This situation is to 
be contrasted with that of the single-particle Green's function, for which one 
knows how to express the expectation values of U and H [21], 02]. Apparently, 
some information has gotten lost in passing from G to 7. This can also be seen 
very clearly from Eq. (|47|) . which shows that information on the dynamics 
of the system, which is contained in G, is erased in the definition of 7(r, r'). 
Explicit information on the static properties of the system is contained in 
the iV-particle density matrix, but as seen from fj4"5l) and (|46j) . a large part of 
this information is also lost ('integrated out') in passing from the iV-particle 
density matrix to the reduced two- or one-particle density matrices. 

Apparently even less information is contained in the particle density^! 
n(r), which is obtained by summing the diagonal element of 7(2, x') over the 
spin variable, 

n{v) = 7(rs, rs). (52) 

s 

This equation follows immediately from comparing (jSJ) with (|46p . We can 
define an alternative density operator, h, by requiring that the same equation 
must also be obtained by substituting n(r) into Eq. (l5Tj) . which holds for any 
single-particle operator. This requirement implies that n(r) = S(r — r,)l 24 l 



22 For nonmultiplicative single-particle operators (such as the kinetic energy, which con- 
tains a derivative) one requires the full single-particle matrix j(x, x') and not only j(x, x). 

23 A quantitative estimate of how much less information is apparently contained in the 
density than in the wave function is given in footnote 6. 

24 The expectation value of n is the particle density and therefore h is often also called 
the density operator. This concept must not be confused with any of the various density 
matrices or the density operator of statistical physics, Eq. |[1D]). 
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The particle density is an even simpler function than j(x,x'): it de- 
pends on one set of coordinates x only, it can easily be visualized as a three- 
dimensional charge distribution, and it is directly accessible in experiments. 
These advantages, however, seem to be more than compensated by the fact 
that one has integrated out an enormous amount of specific information 
about the system in going from wave functions to Green's functions, and on 
to density matrices and finally the density itself. This process is illustrated 
in Fig. HJ 

The great surprise of density-functional theory is that in fact no informa- 
tion has been lost at all, at least as long as one considers the system only 
in its ground state: according to the Hohenberg-Kohn theorem the ground- 
state density n (x) completely determines the ground-state wave function 



Hence, in the ground state, a function of one variable is equivalent to a 
function of N variables! This property shows that we have only integrated out 
explicit information on our way from wave functions via Green's functions 
and density matrices to densities. Implicitly all the information that was 
contained in the ground-state wave function is still contained in the ground- 
state density. Part of the art of practical DFT is how to get this implicit 
information out, once one has obtained the density! 

4 DFT as an effective single-body theory: The 
Kohn-Sham equations 

Density-functional theory can be implemented in many ways. The min- 
imization of an explicit energy functional, discussed up to this point, is not 
normally the most efficient among them. Much more widely used is the Kohn- 
Sham approach. Interestingly, this approach owes its success and popularity 
partly to the fact that it does not exclusively work in terms of the particle (or 
charge) density, but brings a special kind of wave functions (single-particle 
orbitals) back into the game. As a consequence DFT then looks formally 
like a single-particle theory, although many-body effects are still included 

25 The Runge-Gross theorem, which forms the basis of time-dependent DFT [5T], sim- 
ilarly guarantees that the time-dependent density contains the same information as the 
time-dependent wave function. 
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Figure 1: Information on the time-and-space dependent wave func- 
tion ty(x\,X2 ... ,XN,t) is built into Green's functions, and on the time- 
independent wave function into density matrices. Integrating out degrees 
of freedom reduces the iV-particle Green's function and iV-particle density 
matrix to the one-particle quantities G(xi,x 2 ;t) and 7(21, x 2 ) described in 
the main text. The diagonal element of the one-particle density matrix is 
the ordinary charge density — the central quantity in DFT. The Hohenberg- 
Kohn theorem and its time-dependent generalization (the Runge-Gross the- 
orem) guarantee that the densities contain exactly the same information as 
the wave functions. 



via the so-called exchange-correlation functional. We will now see in some 
detail how this is done. 



4.1 Exchange-correlation energy: definition, interpre- 
tation and exact properties 

4.1.1 Exchange-correlation energy 

The Thomas-Fermi approximation ( 1341) for T[n] is not very good. A more 
accurate scheme for treating the kinetic-energy functional of interacting elec- 
trons, T[n], is based on decomposing it into one part that represents the ki- 
netic energy of noninteracting particles of density n, i.e., the quantity called 
above T s [n], and one that represents the remainder, denoted T c [n] (the sub- 
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scripts s and c stand for 'single-particle' and 'correlation', respectively) r 6 l 

T[n]=T s [n]+T c [n}. (53) 

T s [n] is not known exactly as a functional of n [and using the LDA to ap- 
proximate it leads one back to the Thomas- Fermi approximation f)34p ]. but 
it is easily expressed in terms of the single-particle orbitals 0i(r) of a nonin- 
teracting system with density n, as 

h 2 



h 2 N r 
2m ; J 



(54) 



because for noninteracting particles the total kinetic energy is just the sum 
of the individual kinetic energies. Since all fa(r) are functionals of n, this 
expression for T s is an explicit orbital functional but an implicit density 
functional, T s [n] = T s [{fa[n]}], where the notation indicates that T s depends 
on the full set of occupied orbitals fa, each of which is a functional of n. 
Other such orbital functionals will be discussed in Sec. El 
We now rewrite the exact energy functional as 



E[n] = T 


n 


+ U[n] + V 


n 


= T s [{fa[n]}] + U H 


n 


+ E xc 


n] + V[n], 



(55) 



where by definition E xc contains the differences T — T s (i.e. T c ) and U — Uh- 
This definition shows that a significant part of the correlation energy E c 
is due to the difference T c between the noninteracting and the interacting 
kinetic energies. Unlike Eq. (!35"!) . Eq. (|55l) is formally exact, but of course 
E xc is unknown — although the HK theorem guarantees that it is a density 
functional. This functional, E xc [n], is called the exchange-correlation (xc) 
energy. It is often decomposed as E xc = E x + E c , where E x is due to the 
Pauli principle (exchange energy) and E c is due to correlations. (T c is then 
a part of E c .) The exchange energy can be written explicitly in terms of the 
single-particle orbitals aa 27 l 



E x i{fain]}] = - q 4Z 3k fd 3 rfd :i r 



, ^(r)^(r')^(r')«jb(r) 



(56) 



26 T s is defined as the expectation value of the kinetic-energy operator T with the 
Slater determinant arising from density n, i.e., T s [n] = |T|$[n]). Similarly, the full 
kinetic energy is defined as T[n] = (^f[n]\T\^f[n]). All consequences of antisymmetrization 
(i.e., exchange) are described by employing a detcrminantal wave function in defining T s . 
Hence, T C) the difference between T s and T is a pure correlation effect. 

27 This differs from the exchange energy used in Hartree-Fock theory only in the substi- 
tution of Hartree-Fock orbitals <j)f F (r) by Kohn-Sham orbitals <j>i(r). 
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which is known as the Fock term, but no general exact expression in terms 
of the density is known. 



4.1.2 Different perspectives on the correlation energy 

For the correlation energy no general explicit expression is known, neither in 
terms of orbitals nor densities. Different ways to understand correlations are 
described below. 

Correlation energy: variational approach. A simple way to understand 
the origin of correlation is to recall that the Hartree energy is obtained in 
a variational calculation in which the many-body wave function is approx- 
imated as a product of single-particle orbitals. Use of an antisymmetrized 
product (a Slater determinant) produces the Hartree and the exchange energy 
[151 HH] • The correlation energy is then defined as the difference between the 
full ground-state energy (obtained with the correct many-body wave func- 
tion) and the one obtained from the (Hartree-Fock or Kohn-ShanQ) Slater 
determinant. Since it arises from a more general trial wave function than 
a single Slater determinant, correlation cannot raise the total energy, and 
E c [n] < 0. Since a Slater determinant is itself more general than a simple 
product we also have E x < 0, and thus the upper bouncQ E xc [n] < 0. 

Correlation energy: probabilistic approach. Recalling the quantum me- 
chanical interpretation of the wave function as a probability amplitude, we 
see that a product form of the many-body wave function corresponds to treat- 
ing the probability amplitude of the many-electron system as a product of 
the probability amplitudes of individual electrons (the orbitals). Mathemat- 
ically, the probability of a composed event is only equal to the probability 
of the individual events if the individual events are independent (i.e., uncor- 
rected). Physically, this means that the electrons described by the product 
wave function are independent @ Such wave functions thus neglect the fact 
that, as a consequence of the Coulomb interaction, the electrons try to avoid 



28 The Hartree-Fock and the Kohn-Sham Slater determinants are not identical, since 
they are composed of different single-particle orbitals, and thus the definition of exchange 
and correlation energy in DFT and in conventional quantum chemistry is slightly different 
[52]. 

29 A lower bound is provided by the Lieb-Oxford formula, given in Eq. 

30 Correlation is a mathematical concept describing the fact that certain events are not 
independent. It can be defined also in classical physics, and in applications of statistics 
to other fields than physics. Exchange is due to the indistinguishability of particles, and 
a true quantum phenomenon, without any analogue in classical physics. 
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each other. The correlation energy is simply the additional energy lowering 
obtained in a real system due to the mutual avoidance of the interacting 
electrons. One way to characterize a strongly correlated system is to define 
correlations as strong when E c is comparable in magnitude to, or larger than, 
other energy contributions, such as Eh or T s . (In weakly correlated systems 
E c normally is several orders of magnitude smaller .f^l 

Correlation energy: beyond mean-field approach. A rather different (but 
equivalent) way to understand correlation is to consider the following alter- 
native form of the operator representing the Coulomb interaction, equivalent 
to Eq. ©, 

U= q -[d 3 r[ dV Hr)n(r>) -n(r)5(r-r^ 



in which the operator character is carried only by the density operators n 
(in occupation number representation), and the term with the delta function 
subtracts out the interaction of a charge with itself (cf., e.g., the appendix 
of Ref. [53J for a simple derivation of this form of U). The expectation 
value of this operator, U = ($\U\^), involves the expectation value of a 
product of density operators, (\]/|n(r)n(r')|\I>). In the Hartree term (|32|) . 
on the other hand, this expectation value of a product is replaced by a 
product of expectation values, each of the form n(r) = (^^(r)^). This 
replacement amounts to a mean-field approximation, which neglects quantum 
fluctuations^] about the expectation values: by writing h = n + 5hji uc and 
substituting in Eq. (|57|) we see that the difference between (\l/|i7|\l/) and 
the Hartree term (|32|) is entirely due to the fluctuations 5hfi uc and the self- 
interaction correction to the Hartree term. Quantum fluctuations about the 
expectation value are thus at the origin of quantum correlations between 
interacting particles. 



31 Other characterizations of strongly correlated systems are to compare the width of the 
conduction band in a solid with the kinetic energy (if the band width is smaller, correlations 
are strong), or the quasiparticle energies 1i with the Kohn-Sham eigenvalues (if both 
are similar, correlations are weak, see footnote 37), or the derivative discontinuity A xc , 
defined in Eq. (|65j) . with the Kohn-Sham energy gap (if the former is comparable to or 
larger than the latter, correlations are strong). (The meaning of e^, ti and A xc is explained 
below.) No universally applicable definition of 'strong correlations' seems to exist. 

32 At finite temperature there are also thermal fluctuations. To properly include these 
one must use a finite-temperature formulation of DFT [SI] . See also the contribution of 
B. L. Gyorffy et al. in Ref. for DFT treatment of various types of fluctuations. 
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Correlation energy: holes. The fact that both exchange and correlation 
tend to keep electrons apart, has given rise to the concept of an xc hole, 
n xc (r,r'), describing the reduction of probability for encountering an elec- 
tron at r', given one at r. The xc energy can be written as a Hartree-like 
interaction between the charge distribution n(r) and the xc hole n xc (r, r') = 
n x (r,r') + n c (r,r'), 

^M^^/^/dv » (r >;-<;| r ') , (58) 

which defines n xc . The exchange component E x [n] of the exact exchange- 
correlation functional describes the energy lowering due to antisymmetriza- 
tion (i.e., the tendency of like-spin electrons to avoid each other). It gives rise 
to the exchange hole n x (r, r'), which obeys the sum rule / d 3 r' n x (r, r') = — 1. 
The correlation component E c [n] accounts for the additional energy lower- 
ing arising because electrons with opposite spins also avoid each other. The 
resulting correlation hole integrates to zero, so that the total xc hole satis- 
fies / d 3 r' n xc (r, r') = —1. The xc hole can also be written as n xc (r,r') = 
n(r')(g[n](r, r') — 1), where g is the average of the pair-correlation function 
g(r, r'), mentioned in Sec. l3.5.2| over all values of the particle-particle interac- 
tion, from zero (KS system) to (U) (interacting system). This average is sim- 
ply expressed in terms of the coupling constant a as g(r, r') = f g a (r, r')da. 
For the Coulomb interaction, a = e 2 , i.e., the square of the electron charge 

0E]. 

4.1.3 Exact properties 

Clearly E c is an enormously complex object, and DFT would be of little use if 
one had to know it exactly for making calculations. The practical advantage 
of writing E[n] in the form Eq. (|55[) is that the unknown functional -E xc [n] is 
typically much smaller than the known terms T s , Uh and V. One can thus 
hope that reasonably simple approximations for ^^[n] provide useful results 
for E[n}. Some successful approximations are discussed in Sec. [51 Exact prop- 
erties, such as the sum rule / d 3 r'n xc (r, r') = —1, described in the preceding 
section, are most valuable guides in the construction of approximations to 
E xc [n). 

Among the known properties of this functional are the coordinate scaling 
conditions first obtained by Levy and Perdew [55] 

E x [n x ) = XE x [n] (59) 
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E c [n x ] > \E c [n] 
E c [n x ] < XE c [n] 



for A > 1 
for A < 1, 



(60) 
(61) 



where n\(r) = A 3 n(Ar) is a scaled density integrating to total particle number 
N. 

Another important property of the exact functional is the one-electron 
limit 



E c [n {l) ] 
E x [nU] 







-E H [n 



(i)i 



(62) 
(63) 



where vS 1 ' is a one-electron density. These latter two conditions, which are 
satisfied within the Hartree-Fock approximation, but not by standard local- 
density and gradient-dependent functionals, ensure that there is no artificial 
self-interaction of one electron with itself. 
The Lieb-Oxford bound 157], 



E x [n] > E xc [n] > -1.68e 2 J d 3 rn(r) 4/s , 



(64) 



establishes a lower bound on the xc energy, and is satisfied by LDA and many 
(but not all) GGAs. 

One of the most intriguing properties of the exact functional, which has 
resisted all attempts of describing it in local or semilocal approximations, 
is the derivative discontinuity of the xc functional with respect to the total 
particle number 



5E xc [n] 



5n(r) 



SE xc [n] 



N+S 



5n(r) 



(65) 



N-S 



where 5 is an infinitesimal shift of the electron number N, and A xc is a 
system-dependent, but r-independent shift of the xc potential v xc (r) as it 
passes from the electron-poor to the electron-rich side of integer N. The 
noninteracting kinetic-energy functional has a similar discontinuity, given by 



STM 



8n(r) 



N+5 



8n(r) 



e 7V+l — £N 



A 



KS, 



(66) 



N-S 



where ejy and e7v +1 are the Kohn-Sham (KS) single-particle energies of the 
highest occupied and lowest unoccupied eigenstate. The meaning of these 
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KS eigenvalues is discussed in the paragraphs following Eq. (1751) and illus- 
trated in Fig. [2J In the chemistry literature these are called the HOMO 
(highest occupied molecular orbital) and LUMO (lowest unoccupied molec- 
ular orbital), respectively. The kinetic-energy discontinuity is thus sim- 
ply the KS single-particle gap Aks, or HOMO-LUMO gap, whereas the 
xc discontinuity A xc is a many-body effect. The true fundamental gap 
A = E(N + 1) + E(N - 1) - 2E(N) is the discontinuity of the total ground- 
state energy functional [51 [50j E3 [59|, 



A 5E 



<5n(r) 



8 El 



N+S Sn ( r ) 



= A KS + A 3 

N~8 



(67) 



Since all terms in E other than E xc and T s are continuous functionals of 
n(r), the fundamental gap is the sum of the KS gap and the xc discontinuity. 
Standard density functionals (LDA and GGA) predict A xc = 0, and thus 
often underestimate the fundamental gap. The fundamental and KS gaps 
are also illustrated in Fig. [2j 

All these properties serve as constraints or guides in the construction of 
approximations for the functionals E x [n] and £7 c [n]. Many other similar prop- 
erties are known. A useful overview of scaling properties is the contribution 
of M. Levy in Ref. p2]. 



4.2 Kohn-Sham equations 

4.2.1 Derivation of the Kohn-Sham equations 

Since T s is now written as an orbital functional one cannot directly minimize 
Eq. (155]) with respect to n. Instead, one commonly employs a scheme sug- 
gested by Kohn and Sham [BU] for performing the minimization indirectly. 
This scheme starts by writing the minimization as 

Q = SEW = 5T s [n] | 5V[n] | 5U H [n] | 5E xc [n] = 5T s [n] | ^ | 

Sn(r) 5n(r) Sn(r) Sn(r) 5n(r) 5n(r) 

(68) 

As a consequence of Eq. ( l2Tj) . 5V/5n = v(r), the 'external' potential the 
electrons move injff] The term SUn/Sn simply yields the Hartree potential, 



33 This potential is called 'external' because it is external to the electron system and 
not generated self-consistently from the electron-electron interaction, as vh and v xc . It 
comprises the lattice potential and any additional truly external field applied to the system 
as a whole. 
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introduced in Eq. fill I) . For the term 5E xc /6n, which can only be calculated 
explicitly once an approximation for E xc has been chosen, one commonly 
writes v xc . By means of the Sham-Schliiter equation fl39|) . v xc is related to 
the irreducible self energy £, introduced in Eq. fl38l) [50J. 

Consider now a system of noninteracting particles moving in the potential 
v s (r). For this system the minimization condition is simply 







SE s [n] _ 5T s [n] 5V s [n] _ 6T s [n 



+ 



5n(r) 5n(r) 5n(r) 8n(r 



+ v 3 (r) 



(69) 



since there are no Hartree and xc terms in the absence of interactions. The 
density solving this Euler equation is n s (r). Comparing this with Eq. (1681) 
we find that both minimizations have the same solution n s (r) = n(r), if v s 
is chosen to be 

' ' (70) 



v a (r) = v(r) + v H (r) + v xc (r). 



Consequently, one can calculate the density of the interacting (many-body) 
system in potential v(r), described by a many-body Schrodinger equation 
of the form (J2]), by solving the equations of a noninteracting (single-body) 
system in potential u s (r)o 

In particular, the Schrodinger equation of this auxiliary system, 



ft 2 V 2 
2m 



(71) 



yields orbitals that reproduce the density n(r) of the original system (these 
are the same orbitals employed in Eq. fET 



(72) 



n r 



n s r 



where /j is the occupation of the i'th orbitalQ Eqs. (I7D|) to (1721) are the 
celebrated Kohn-Sham (KS) equations. They replace the problem of mini- 
mizing E[n] by that of solving a noninteracting Schrodinger equation. (Recall 



34 The question whether such a potential v s (r) always exists in the mathematical sense 
is called the noninteracting w-representability problem. It is known that every interacting 
ensemble u-representable density is also noninteracting ensemble i;-representable, but, as 
mentioned in Sec. 13.21 only in discretized systems has it been proven that all densities are 
interacting ensemble w-representable. It is not known if interacting ensemble-representable 
densities may be noninteracting pure-state representable (i.e, by a single determinant), 
which would be convenient (but is not necessary) for Kohn-Sham calculations. 

35 Normally, the occupation numbers /j follows an Aufbau principle (Fermi statistics) 
with fa = 1 for i < N, fa = for i > N, and < fa < 1 for i — N (i.e., at most 
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that the minimization of E[n] originally replaced the problem of solving the 
many-body Schrodinger equation!) 

Since both vh and v xc depend on n, which depends on the (pi, which in 
turn depend on v s , the problem of solving the KS equations is a nonlinear 
one, just as is the one of solving the (much more complicated) Dyson equation 
(|38p . The usual way of solving such problems is to start with an initial guess 
for n(r), calculate the corresponding v s (r), and then solve the differential 
equation (ITT!) for the 4>i- From these one calculates a new density, using (!72|) . 
and starts again. The process is repeated until it converges. The technical 
name for this procedure is 'self-consistency cycle'. Different convergence 
criteria (such as convergence in the energy, the density, or some observable 
calculated from these) and various convergence-accelerating algorithms (such 
as mixing of old and new effective potentials) are in common use. Only rarely 
it requires more than a few dozen iterations to achieve convergence, and even 
rarer are cases where convergence seems unattainable, i.e., a self-consistent 
solution of the KS equation cannot be found. 

Once one has a converged solution n , one can calculate the total energy 
from Eq. f[53|) or, equivalent ly and more conveniently, fronj^l 



E = Ef U - % I rf 3 r / dV - J d 3 r ^ c (r)n (r) + E xc [n ). 



(73) 



Equation f!73|) follows from writing V[n] in fl55l) by means of flTOl) as 



V\n] 



J d 3 rt)(r)n(r) = J d 3 r [v s (r) — w#(r) — v xc (r)]n(r) (74) 
= V,[n] - / d 3 r [v H (r) + v xc (r)]n(r), (75) 



the highest occupied orbital can have fractional occupation). Some densities that are not 
nonintcracting w-representable by a single ground-state Slater determinant, may still be 
obtained from a single determinant if one uses occupation numbers f L that leave holes 
below the HOMO (the Fermi energy in a metal), so that fi ^ 1 even for some i < N 
[5T] . but this is not guaranteed to describe all possible densities. Alternatively (see Sec. 
3.21 and footnote 34) a Kohn-Sham equation may be set up in terms of ensembles of 
determinants. This guarantees noninteracting u-representability for all densities that are 
interacting ensemble u-representable. For practical KS calculations, the most important 
consequence of the fact that not every arbitrary density is guaranteed to be noninteracting 
w-representable is that the Kohn-Sham sclfconsistency cycle is not guaranteed to converge. 

36 All terms on the right-hand side of (f73|) except for the first, involving the sum of the 
single-particle energies, are sometimes known as double-counting corrections, in analogy 
to a similar equation valid within Hartree-Fock theory. 
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and identifying the energy of the noninteracting (Kohn-Sham) system as 
E. = Y?*i = T a + V a . 



4.2.2 The eigenvalues of the Kohn-Sham equation 

Equation f!73|) shows that E is not simply the surrj^l of all e^. In fact, it 
should be clear from our derivation of Eq. (ITT]) that the e; are introduced as 
completely artificial objects: they are the eigenvalues of an auxiliary single- 
body equation whose eigenfunctions (orbitals) yield the correct density It is 
only this density that has strict physical meaning in the KS equations. The 
KS eigenvalues, on the other hand, in general bear only a semiquantitative 
resemblance with the true energy spectrum [61], but are not to be trusted 
quantitatively 

The main exception to this rule is the highest occupied KS eigenvalue. 
Denoting by €n(M) the iV'th eigenvalue of a system with M electrons, one 
can show rigorously that €n{N) = —I, the negative of the first ionization 
energy of the iV-body system, and e7v + i(iV + 1) = —A, the negative of the 
electron affinity of the same iV-body system [58], [63]. These relations 
hold for the exact functional only. When calculated with an approximate 
functional of the LDA or GGA type, the highest eigenvalues usually do not 
provide good approximations to the experimental I and A. Better results 
for these observables are obtained by calculating them as total-energy differ- 
ences, according to I = E (N - 1) - E (N) and A = E (N) — E (N + 1), 
where E (N) is the ground-state energy of the iV-body system. Alterna- 
tively, self-interaction corrections can be used to obtain improved ionization 
energies and electron affinities from Kohn-Sham eigenvalues [64|. 

Figure [2] illustrates the role played by the highest occupied and lowest 
unoccupied KS eigenvalues, and their relation to observables. For molecules, 
HOMO(N) is the highest-occupied molecular orbital of the iV-electron sys- 
tem, HOMO(N+l) that of the N + 1-electron system, and LUMO(N) the 
lowest unoccupied orbital of the iV-electron system. In solids with a gap, 

37 The difference between Eq and e, is due to particle-particle interactions. The ad- 
ditional terms on the right-hand side of (|T3[) give mathematical meaning to the common 
statement that the whole is more than the sum of its parts. If Eq can be written approx- 
imately as Ylf <U (where the <f, are not the same as the KS eigenvalues ej) the system 
can be described in terms of N weakly interacting quasiparticles, each with energy e^. 
Fermi-liquid theory in metals and effective-mass theory in semiconductors are examples 
of this type of approach. 
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Figure 2: Schematic description of some important Kohn-Sham eigenvalues 
relative to the vacuum level, denoted by 0, and their relation to observables. 
See main text for explanations. 

the HOMO and LUMO become the top of the valence band and the bottom 
of the conduction band, respectively, whereas in metals they are both iden- 
tical to the Fermi level. The vertical lines indicate the Kohn-Sham (single- 
particle) gap Aks, the fundamental (many-body) gap A, the derivative dis- 
continuity of the xc functional, A xc , the ionization energy of the interacting 
iV-electron system I(N) = — e^^N) (which is also the ionization energy of 
the Kohn-Sham system Iks{N)), the electron affinity of the interacting iV- 
electron system A(N) = —eN + i(N + 1) and the Kohn-Sham electron affinity 
Aks(N) = -e N+1 (N). 

Given the auxiliary nature of the other Kohn-Sham eigenvalues, it comes 
as a great (and welcome) surprise that in many situations (typically char- 
acterized by the presence of fermionic quasiparticles and absence of strong 
correlations) the Kohn-Sham eigenvalues e» do, empirically, provide a rea- 
sonable first approximation to the actual energy levels of extended systems. 
This approximation is behind most band-structure calculations in solid-state 
physics, and often gives results that agree well with experimental photoemis- 
sion and inverse photoemission data [65], but much research remains to be 
done before it is clear to what extent such conclusions can be generalized, 
and how situations in which the KS eigenvalues are good starting points for 
approximating the true excitation spectrum are to be characterized micro- 
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scopically PJETjB 

Most band-structure calculations in solid-state physics are actually calcu- 
lations of the KS eigenvalues e,Jf§ This simplification has proved enormously 
successful, but when one uses it one must be aware of the fact that one is 
taking the auxiliary single-body equation (ITT]) literally as an approximation 
to the many-body Schrodinger equation. DFT, practiced in this mode, is not 
a rigorous many-body theory anymore, but a mean-field theory (albeit one 
with a very sophisticated mean field v s (r)). 

The energy gap obtained in such band-structure calculations is the one 
called HOMO-LUMO gap in molecular calculations, i.e., the difference be- 
tween the energies of the highest occupied and the lowest unoccupied single- 
particle states. Neglect of the derivative discontinuity A xc , defined in Eq. (1B3D . 
by standard local and semilocal xc functionals leads to an underestimate of 
the gap (the so-called 'band-gap problem'), which is most severe in transition- 
metal oxides and other strongly correlated systems. Self-interaction correc- 
tions provide a partial remedy for this problem [7T1 U2i E3 El]- 



4.2.3 Hartree, Hartree-Fock and Dyson equations 



A partial justification for the interpretation of the KS eigenvalues as start- 
ing point for approximations to quasi-particle energies, common in band- 
structure calculations, can be given by comparing the KS equation with 
other self-consistent equations of many-body physics. Among the simplest 
such equations are the Hartree equation 

+ v(r) + v H {r) 



2m 

and the Hartree-Fock (HF) equation 
+ v(r) + v H (r 



kHz 



WW 



(76) 
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2m 



lHF i 



(I 



3./ 7(r, 



d r 



(77) 

where 7(r, r') is the density matrix of Eq. (14*61) . It is a fact known as Koop- 
man's theorem [49J that the HF eigenvalues ef F can be interpreted as unre- 
laxed electron-removal energies (i.e., ionization energies of the z'th electron, 



38 Several more rigorous approaches to excited states in DFT, which do not require the 
KS eigenvalues to have physical meaning, are mentioned in Sec. [6] 

A computationally more expensive, but more reliable, alternative is provided by the 
GW approximation [Ml [69l EH] • 
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neglecting reorganization of the remaining electrons after removal). As men- 
tioned above, in DFT only the highest occupied eigenvalue corresponds to 
an ionization energy, but unlike in HF this energy includes relaxation effects. 

The KS equation fITTT) includes both exchange and correlation via the 
multiplicative operator v xc . Both exchange and correlation are normally 
approximated in DFT0 whereas HF accounts for exchange exactly, through 
the integral operator containing 7(r, r'), but neglects correlation completely. 
In practise DFT results are typically at least as good as HF ones and often 
comparable to much more sophisticated correlated methods — and the KS 
equations are much easier to solve than the HF equations]^] 

All three single-particle equations, Hartree, Hartree-Fock and Kohn-Sham 
can also be interpreted as approximations to Dyson's equation (1381) . which 
can be rewritten as [H] 

2V72 



2m 



+ v(r) j V*(r) + J d 3 r' £(r, r', E k )^ k {r') = E k ^ k {r), (78) 



where £ is the irreducible self energy introduced in Eq. fl38|) . The E k ap- 
pearing in this equation are the true (quasi-)electron addition and removal 
energies of the many-body system. Needless to say, it is much more compli- 
cated to solve this equation than the HF or KS equations. It is also much 
harder to find useful approximations for E than for f xc @ Obviously, the KS 
equation employs a local, energy-independent potential v s in place of the non- 
local, energy- dependent operator S. Whenever this is a good approximation, 
the ei are also a good approximation to the E k . 

The interpretation of the KS equation ( JTB as an approximation to Eq. (17H1) 
is suggestive and useful, but certainly not necessary for DFT to work: if the 
KS equations are only used to obtain the density, and all other observables, 
such as total energies, are calculated from this density, then the KS equa- 
tions in themselves are not an approximation at all, but simply a very useful 
mathematical tool. 



40 A possibility to treat exchange exactly in DFT is offered by the OEP method discussed 
in Sec. 15.31 

41 This is due to the integral operator in the HF equations. 

42 The GW approximation 68, 69j[70], mentioned in footnote 39, is one such approxi- 
mation for S, but in actual implementations of it one usually takes DFT-KS results as an 
input. 
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4.3 Basis functions 

In practice, numerical solution of the KS differential equation flTTj) typically 
proceeds by expanding the KS orbitals in a suitable set of basis functions and 
solving the resulting secular equation for the coefficients in this expansion 
and/or for the eigenvalues for which it has a solution. The construction 
of suitable basis functions is a major enterprise within electronic-structure 
theory (with relevance far beyond DFT), and the following lines do little 
more than explaining some acronyms often used in this field. 

In physics much is known about the construction of basis functions for 
solids due to decades of experience with band-structure calculations. This in- 
cludes many calculations that predate the widespread use of DFT in physics. 
There is a fundamental dichotomy between methods that work with fixed ba- 
sis functions that do not depend on energy, and methods that employ energy- 
dependent basis functions. Fixed basis functions are used e.g., in plane- wave 
expansions, tight-binding or LCAO (linear combination of atomic orbitals) 
approximations, or the OPW (orthogonalized plane wave) method. Exam- 
ples for methods using energy- dependent functions are the AP W (augmented 
plane wave) or KKR (Korringa-Kohn-Rostoker) approaches. This distinction 
became less clear-cut with the introduction of 'linear methods' [75], in which 
energy-dependent basis functions are linearized (Taylor expanded) around 
some fixed reference energy. The most widely used methods for solving the 
Kohn-Sham equation in solid-state physics, LMTO (linear muffin tin or- 
bitals) and LAPW (linear augmented plane waves), are of this latter type 
[76] . Development of better basis functions is an ongoing enterprise (TTJ [78]. 

The situation is quite similar in chemistry. Due to decades of experience 
with Hartree-Fock and CI calculations much is known about the construction 
of basis functions that are suitable for molecules. Almost all of this continues 
to hold in DFT — a fact that has greatly contributed to the recent popularity 
of DFT in chemistry. Chemical basis functions are classified with respect 
to their behaviour as a function of the radial coordinate into Slater type 
orbitals (STOs), which decay exponentially far from the origin, and Gaussian 
type orbitals (GTOs), which have a gaussian behaviour. STOs more closely 
resemble the true behaviour of atomic wave functions [in particular the cusp 
condition of Eq. (fl9~j) ]. but GTOs are easier to handle numerically because 
the product of two GTOs located at different atoms is another GTO located 
in between, whereas the product of two STOs is not an STO. The so-called 
'contracted basis functions', in which STO basis functions are reexpanded in 
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a small number of GTOs, represent a compromise between the accuracy of 
STOs and the convenience of GTOs. The most common methods for solving 
the Kohn-Sham equations in quantum chemistry are of this type [H H9]. 
Very accurate basis functions for chemical purposes have been constructed 
by Dunning [79] and, more recently, by da Silva and collaborators [80, [81] . 
More details on the development of suitable basis functions can be found, 
e.g., in these references and Ref. |49j . 

A very popular approach to larger systems in DFT, in particular solids, 
is based on the concept of a pseudopotential (PP). The idea behind the PP 
is that chemical binding in molecules and solids is dominated by the outer 
(valence) electrons of each atom. The inner (core) electrons retain, to a good 
approximation, an atomic-like configuration, and their orbitals do not change 
much if the atom is put in a different environment. Hence, it is possible to 
approximately account for the core electrons in a solid or a large molecule 
by means of an atomic calculation, leaving only the valence density to be 
determined self-consistently for the system of interest. 

In the original Kohn-Sham equation the effective potential v s [n] = v ext + 
vhV 1 ] + UsscM is determined by the full electronic density n(r), and the self- 
consistent solutions are single-particle orbitals reproducing this density. In 
the PP approach the Hartree and xc terms in v s [n] are evaluated only for 
the valence density n v , and the core electrons are accounted for by replac- 
ing the external potential v ext by a pseudopotential v pp t . Hence v pp [n v ] = 
v eJt + v H[ n v] + Va;c[w«]o The PP v p ^ is determined in two steps. First, one 
determines, in an auxiliary atomic calculation, an effective PP, v pp , such 
that for a suitably chosen atomic reference configuration the single-particle 
orbitals resulting from v P p agree — outside a cut-off radius r c separating the 
core from the valence region — with the valence orbitals obtained from the 
all-electron KS equation for the same atom. As a consequence, the valence 
densities n"* obtained from the atomic KS and the atomic PP equation are 
the same. Next, one subtracts the atomic valence contributions f//[n°*] and 
Vxcin^f] from v p p [n% t ] to obtain the external PP f^f J3 which is then used 

43 Note that the effective potential v s is a way to deal with the electron-electron inter- 
action. The pseudopotential is a way to deal with the density of the core electrons. Both 
potentials can be profitably used together, but are conceptually different. 

44 This external PP is also called the unscreened PP, and the subtraction of Ufl-fn"*] and 
v xc[ n v t ] from «f is called the 'unscreening of the atomic PP'. It can only be done 

exactly for the Hartree term, because the contributions of valence and core densities are 
not additive in the xc potential (which is a nonlinear functional of the total density) . 
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in the molecular or solid-state calculation, together with %[nj and u^n,,] 
taken at the proper valence densities for these systems. 

The way v pp is generated from the atomic calculation is not unique. 
Common pseudopotentials are generated following the prescription of, e.g., 
Bachelet, Hamann and Schliiter [82], Kleinman and Bylander [83J, Vanderbilt 
[HI] or Troullier and Martins [BS]. Useful reviews are Refs. [SSI E3 EH] ■ The 
pseudopotential approach is very convenient because it reduces the number of 
electrons treated explicitly, making it possible to perform density-functional 
calculations on systems with tens of thousands of electrons. Moreover, the 
pseudopotentials v p x P are much smoother than the bare nuclear potentials 
v ext . The remaining valence electrons are thus well described by plane- wave 
basis sets. 

Some of the choices one has to make in a practical Kohn-Sham calculation 
are illustrated schematically in Fig. [3j 

5 Making DFT practical: Approximations 

There are basically three distinct types of approximations involved in 
a DFT calculation. One is conceptual, and concerns the interpretation of 
KS eigenvalues and orbitals as physical energies and wave functions. This 
approximation is optional — if one does not want to make it one simply 
does not attach meaning to the eigenvalues of Eq. (ITTj) . The pros and cons 
of this procedure were discussed in Sees. 14.2.21 and 14.2.31 The second type 
of approximation is numerical, and concerns methods for actually solving 
the differential equation (1711) . A main aspect here is the selection of suitable 
basis functions, briefly discussed in Sec. 14.31 The third type of approximation 
involves constructing an expression for the unknown xc functional i? xc [n], 
which contains all many-body aspects of the problem [cf. Eq. (1551) ]. It is 
with this type of approximation that we are concerned in the present section. 

This chapter is intended to give the reader an idea of what types of 
functionals exist, and to describe what their main features are, separately 
for local functionals (TF, LDA and Xa; Sec. 15.11) . semilocal, or gradient- 
dependent, functionals (GEA and GGA; Sec. I5.2[) . and nonlocal functionals 
(hybrids, orbital functionals such as meta-GGAs, EXX and SIC, and integral- 
dependent functionals such as ADA; Sec. 15.31) . This chapter does deal only 
most superficially with the actual construction of these functionals. For more 
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Figure 3: Some of the choices made in a Kohn-Sham calculation. The treat- 
ment can be nonrelativistic (based on Schrodinger's equation), scalar rela- 
tivistic (using the relativistic kinetic-energy operator and other simple rela- 
tivistic corrections, but no spin-orbit coupling) or relativistic (using Dirac's 
equation, which includes also spin-orbit coupling). The core electrons can 
be treated explicitly (all electron calculation) or incorporated, together with 
v ext , in a pseudopotential. The Hartree potential can be obtained from in- 
tegrating the charge density or from solving Poisson's differential equation. 
Many choices are available for the xc potential. The eigenvalues can be 
considered mere Lagrange multipliers or interpreted as zero-order approxi- 
mations to the actual energy spectrum. The eigenfunctions can similarly be 
considered auxiliary functions generating the density, or interpreted as zero- 
order approximations to quasi-particle wave functions. Solution of the KS 
equation can proceed on a numerical mesh, or by expansion of the eigenfunc- 
tions in basis functions. Many types of suitable basis functions exist. For 
every new problem a suitable combination of choices must be made, and all 
possibilities continue to be useful and to be actively explored in physics and 
chemistry. 
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details on functional construction and testing the reader is referred to the 
reviews [5-19] or to the original papers cited below. Sticking to the bird's-eye 
philosophy of this overview I have also refrained from including numerical 
data on the performance of each functional — extensive comparisons of a 
wide variety of functionals can be found in Refs. [5-19] and in the original 
literature cited below. 

5.1 Local functionals: LDA 

Historically (and in many applications also practically) the most important 
type of approximation is the local-density approximation (LDA). To under- 
stand the concept of an LDA recall first how the noninter acting kinetic en- 
ergy T s [n] is treated in the Thomas- Fermi approximation: In a homogeneous 
system one knows that, per volumj^l 

3fi 2 

thorny = ^(3^)2/3^5/3 (79) 

10m 

where n = const. In an inhomogeneous system, with n = n(r), one approxi- 
mates locally 

3fi 2 

t s (v) « t h s ° m (n(r)) = ^( 37 r 2 ) 2 / 3 n(r) 5 / 3 (80) 

10m 

and obtains the full kinetic energy by integration over all space 

T s LDA [n] = J <Prt h ° m {n{v)) = ^-(3tt 2 ) 2 / 3 J d 3 rn(r) 5 / 3 . (81) 

For the kinetic energy the approximation T s [n] m T^ DA [n] is much inferior 
to the exact treatment of T s in terms of orbitals, offered by the Kohn-Sham 
equations, but the LDA concept turned out to be highly useful for another 
component of the total energy fl55l) . the exchange- correlation energy i? xc [n]. 
For the exchange energy E x [n] the procedure is very simple, since the per- 
volume exchange energy of the homogeneous electron liquid is known exactly 

13 El, 

er>) = -^(£) 1/3 n 4 / 3 , (82) 



45 The change from a capital T to a lower-case t is commonly used to indicate quantities 
per volume. 
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so that 

_ _ . „ o , _\ l /a 

(83) 



^ A W = -¥(f) 1/3 /rf 3 m(r) 4 / 3 . 



This is the LDA for £ x O 

For the correlation energy E c [n] the situation is more complicated since 
(n) is not known exactly: the determination of the correlation energy 
of a homogeneous interacting electron system (an electron liquid) is already 
a difficult many-body problem on its own! Early approximate expressions 
for e^ om (n) were based on applying perturbation theory (e.g. the random- 
phase approximation) to this problem [891190] . These approximations became 
outdated with the advent of highly precise Quantum Monte Carlo (QMC) 
calculations for the electron liquid, by Ceperley and Alder [HI]. Modern 
expressions for e^ om (n) [92, 93, 94J are parametrizations of these data. These 
expressions are implemented in most standard DFT program packages and 
in typical applications give almost identical results. On the other hand, the 
earlier parametrizations of the LDA, based on perturbation theory [89, 90J, 
can occasionally deviate substantially from the QMC ones, and are better 
avoided. 

Independently of the parametrization, the LDA for S xc [n] formally con- 
sists iio 



E xc [n] « E™ A [n] = I d 3 re^(n)\ n ^ n{l) = J d 3 re^(n(r)), (84) 



where e^° m = e h ° m + e^ om . The corresponding xc potential is simply 



a horn/ 

vL DA [n}(r) = ^ H 



dn 



(85) 



n— >n(r) 



This approximation for -E xc [n] has proved amazingly successful, even when 
applied to systems that are quite different from the electron liquid that forms 
the reference system for the LDA. A partial explanation for this success of the 



46 If one adds this term to the Thomas-Fermi expression ([35]) one obtains the so-called 
Thomas- Fermi-Dirac approximation to E[n]. It one multiplies it with an adjustable pa- 
rameter a one obtains the so-called Xa approximation to E xc [n]. These approximations 
are not much used today in DFT. 

47 Sometimes one uses the per-particle instead of the per- volume energy of the homoge- 
neous system in writing the LDA. Since the conversion factor between both is the number 
of particles per volume, i.e., the density, an additional n(r) then appears under the inte- 
grals in (fM| and also contributes to 
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LDA is systematic error cancellation: typically, LDA underestimates E c but 
overestimates E x , resulting in unexpectedly good values of E xc . This error 
cancellation is not accidental, but systematic, and caused by the fact that for 
any density the LDA xc hole satisfies the correct sum rule / d 3 r'n x ^ A (r, r') = 
— 1 (see Sec. 14.1.21) . which is only possible if integrated errors in n x DA cancel 
with those of n^ DA . 

For many decades the LDA has been applied in, e.g., calculations of band 
structures and total energies in solid-state physics. In quantum chemistry 
it is much less popular, because it fails to provide results that are accurate 
enough to permit a quantitative discussion of the chemical bond in molecules 
(so-called 'chemical accuracy' requires calculations with an error of not more 
than about 1 kcal/mol = 0.04336 eV/particle). 

At this stage it may be worthwhile to recapitulate what practical DFT 
does, and where the LDA enters its conceptual structure: What real systems, 
such as atoms, molecules, clusters and solids, have in common, is that they 
are simultaneously inhomogeneous (the electrons are exposed to spatially 
varying electric fields produced by the nuclei) and interacting (the electrons 
interact via the Coulomb interaction). The way density- functional theory, in 
the local-density approximation, deals with this inhomogeneous many-body 
problem is by decomposing it into two simpler (but still highly nontrivial) 
problems: the solution of a spatially homogeneous interacting problem (the 
homogeneous electron liquid) yields the uniform xc energy e x ° m (n), and the 
solution of a spatially inhomogeneous noninteracting problem (the inhomoge- 
neous electron gas described by the KS equations) yields the particle density. 
Both steps are connected by the local-density potential fl85|) . which shows how 
the xc energy of the uniform interacting system enters the equations for the 
inhomogeneous noninteracting system. 

The particular way in which the inhomogeneous many-body problem is 
decomposed, and the various possible improvements on the LDA, are behind 
the success of DFT in practical applications of quantum mechanics to real 
materials. Some such improvements on the LDA are discussed in the next 
two sections. 

5.2 Semilocal functionals: GEA, GGA and beyond 

In the LDA one exploits knowledge of the density at point r. Any real system 
is spatially inhomogeneous, i.e., it has a spatially varying density n(r), and 
it would clearly be useful to also include information on the rate of this 
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variation in the functional. A first attempt at doing this were the so-called 
'gradient-expansion approximations' (GEA). In this class of approximation 
one tries to systematically calculate gradient-corrections of the form |Vra(r)|, 
|Vn(r)| 2 , V 2 n(r), etc., to the LDA. A famous example is the lowest-order 
gradient correction to the Thomas- Fermi approximation for T s [n], 

T s [n] ~ Tf [n] = T^[n] + f f d\ l^C. (86) 

8m J n[r) 

This second term on the right-hand side is called the Weizsacker termini 
Similarly, in 

E \n] » E GEAi - 2) \n\ = E LDA \n] — ^ / d 3 r \ Vn ^ (g 7 ) 

x[ 1 x M x [ 1 4327r(37r2)i/3 J a n ( r )4/3 

the second term on the right-hand side is the lowest-order gradient correc- 
tion^ to E^ DA [n\. In practice, the inclusion of low-order gradient corrections 
almost never improves on the LDA, and often even worsens it. Higher-order 
corrections (e.g., oc |Vn(r)| a or oc V /3 n(r) with a,/3 > 2), on the other hand, 
are exceedingly difficult to calculate, and little is known about them. 

In this situation it was a major breakthrough when it was realized, in the 
early eighties, that instead of power-series-like systematic gradient expan- 
sions one could experiment with more general functions of n(r) and Vn(r), 
which need not proceed order by order. Such functionals, of the general form 

Eg GA [n} = /tfV/Hr),Vn(r)), (88) 

have become known as generalized-gradient approximations (GGAs) (95]. 

Different GGAs differ in the choice of the function f(n, Vn). Note that 
this makes different GGAs much more different from each other than the 
different parametrizations of the LDA: essentially there is only one correct 

48 If one adds this term to the Thomas- Fermi expression ([55)) one obtains the so-called 
Thomas- Fermi- Weizsacker approximation to E[n]. In a systematic gradient expansion the 
8 in the denominator is replaced by a 72 [3 [6] . 

49 Remarkably, the form of this term is fully determined already by dimensional analysis: 
In Ex EA ^ = q 2 J d 3 r f(n, |Vri| 2 ) the function / must have dimensions (length)" 4 . Since 
the dimensions of n and |Vn| 2 are (length) -3 and (length)" 8 , respectively, and to second 
order no higher powers or higher derivatives of n are allowed, the only possible combination 
is / oc |Vn(r)| 2 /n 4 / 3 . 
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expression for e^° m (n), and the various parametrizations of the LDA [HHl 
inni E2J E2 El] are merely different ways of writing it. On the other hand, 
depending on the method of construction employed for obtaining f(n, Vn) 
one can obtain very different GGAs. In particular, GGAs used in quantum 
chemistry typically proceed by fitting parameters to test sets of selected 
molecules. On the other hand, GGAs used in physics tend to emphasize 
exact constraints. Nowadays the most popular (and most reliable) GGAs 
are PBE (denoting the functional proposed in 1996 by Perdew, Burke and 
Ernzerhof [96]) in physics, and BLYP (denoting the combination of Becke's 
1988 exchange functional [97J with the 1988 correlation functional of Lee, 
Yang and Parr [HE]) in chemistry. Many other GGA-type functionals are 
also available, and new ones continue to appear. 

Quite generally, current GGAs seem to give reliable results for all main 
types of chemical bonds (covalent, ionic, metallic and hydrogen bridge). For 
van der Waals interactions, however, common GGAs and LDA failO To de- 
scribe these very weak interactions several more specialized approaches have 
been developed within DFT [Ml [1021 EU EDI IDS]. Both in physics and in 
chemistry the widespread use of GGAs has lead to major improvements as 
compared to LDA. 'Chemical accuracy', as defined above, has not yet been 
attained, but is not too far away either. A useful collection of explicit expres- 
sions for some GGAs can be found in the appendix of Ref. [106] . and more 
detailed discussion of some selected GGAs and their performance is given in 
Ref. [lU7j and in the chapter of Kurth and Perdew in Refs. [THl [TB] . 

No systematic attempt at comparing explicit functionals can be made 
here, but many detailed comparisons are available in the literature. For 
pure illustrative purposes only, Table [1] contains ground-state energies of the 
Ar atom, obtained with several of the methods discussed previously in this 
chapter. Footnote 7 contains additional information on the performance of 
DFT for larger systems. 

50 The PBE GGA [H] and the TPSS MGGA g] (see below) may be partial exceptions 
[991 II 00] because they work reasonably well near the equilibrium distance of the van der 
Waals bond, but they recover only the short-range behaviour and do not describe correctly 
the long-range asymptotic regime of the van der Waals interaction. 
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method 


-E/a.u. 


Thomas-Fermi 


625.7 


Hartree-Fock 


526.818 


OEP (exchange only) 


526.812 


LDA (exchange only) 


524.517 


LDA (VWN) 


525.946 


LDA (PW92) 


525.940 


LDA-SIC(PZ) 


528.393 


ADA 


527.322 


WDA 


528.957 


GGA (B88LYP) 


527.551 


experiment 


527.6 



Table 1: Ground-state energy in atomic units (1 a.u. =1 Hartree = 2 
Rydberg = 27.21e^=627.5/cca//mo/) of the Ar atom (Z = 18), obtained with 
some representative density functionals and related methods. The Hartree- 
Fock and OEP (exchange only) values are from Krieger et al. (third of Ref. 
[120J), ADA and WDA values are from Gunnarsson et al., Ref. [129J, as 
reported in Ref. [5j, and the LDA-SIC(PZ) value is from Perdew and Zunger, 
Ref. [93J. The experimental value is based on Veillard and Clementi, J. Chem. 
Phys. 49, 2415 (1968), and given to less significant digits than the calculated 
values, because of relativistic and quantum electrodynamical effects (Lamb 
shift) that are automatically included in the experimental result but not in 
the calculated values. 
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5.3 Orbital functionals and other nonlocal approxima- 
tions: hybrids, Meta-GGA, SIC, OEP, etc. 

In spite of these advances, the quest for more accurate functionals goes ever 
on, and both in chemistry and physics various beyond-GGA functionals have 
appeared. Perhaps the most popular functional in quantum chemistry^! is 
B3LYP. This is a combination of the LYP GGA for correlation [98] with 
Becke's three-parameter hybrid functional B3 for exchange [108J. Common 
hybrid functionals, such as B3, mix a fraction of Hartree-Fock exchange into 
the DFT exchange functional (other mixtures are also possible). The con- 
struction of hybrid functional involves a certain amount of empiricism in the 
choice of functionals that are mixed and in the optimization of the weight 
factors given to the HF and DFT terms. Formally, this might be considered a 
drawback, but in practice B3 has proven to be the most successful exchange 
functional for chemical applications, in particular when combined with the 
LYP GGA functional for E c . More extreme examples of this semiempirical 
mode of construction of functionals are Becke's 1997 hybrid functional [109J, 
which contains 10 adjustable parameters, and the functionals of Refs. jl 10] 
and each of which contains 21 parameters. 

Another recent beyond-GGA development is the emergence of so-called 
Meta-GGAs, which depend, in addition to the density and its derivatives, 
also on the Kohn-Sham kinetic-energy density r(r) [21 11121 1113] 

^) = ^ElV^(r)r, (89) 

so that E xc can be written as £ w [n(r),Vn(r),r(r)]. The additional degree 
of freedom provided by r is used to satisfy additional constraints on E xc , 
such as a self-interaction-corrected correlation functional, recovery of the 
fourth-order gradient expansion for exchange in the limit of slowly varying 
densities, and a finite exchange potential at the nucleus [2]. In several recent 
tests [21 EH EH ESI UTS] Meta-GGAs have given favorable results, even 
when compared to the best GGAs, but the full potential of this type of 
approximation is only beginning to be explored systematically. 

As we have seen in the case of T s , it can be much easier to represent a 
functional in terms of single-particle orbitals than directly in terms of the 

51 This was written in early 2002, but at the time of revision of this text in 2006 it is 
still correct. 
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density. Such functionals are known as orbital functionals, and Eq. (1541) con- 
stitutes a simple example. Another important orbital- dependent functional 
is the exchange energy (Fock term) of Eq. fl56|) . The Meta-GGAs and hybrid 
functionals mentioned above are also orbital functionals, because they de- 
pend on the kinetic energy density fl89l) . and on a combination of the orbital 
functional (15T)|) with ordinary GGAs, respectively. 

Still another type of orbital functional is the self-interaction correction 
(SIC). Most implementations of SIC make use of the expressions proposed 
in Ref. [93] (PZ-SIC), 

E a r ox,SIC [nh ^ = E ap P ro X[nh ^ _ £ ^ ^ _ JjgP".^ (,]) , (90) 

which subtracts, orbital by orbital, the contribution the Hartree and the xc 
functionals would make if there was only one electron in the system. This 
correction can be applied on top of any approximate density functional, and 
ensures that the resulting corrected functional satisfies E^ prox ' SIC [n^ 1 \0] = 
-E H [nW] for a one-electron system. The LDA is exact for a completely 
uniform system, and thus is self- interaction free in this limit, but neither it 
nor common GGAs satisfy the requirement of freedom from self-interaction 
in general, and even Meta-GGAs have a remaining self-interaction error in 
their exchange part [21 1112j . This self- interaction is particularly critical for 
localized states, such as the d states in transition-metal oxides. For such 
systems PZ-SIC has been shown to greatly improve the uncorrected LDA 
[TTj [72] . but for thermochemistry PZ-SIC does not seem to be significant 

Unfortunately the PZ-SIC approach, which minimizes the corrected en- 
ergy functional with respect to the orbitals, does not lead to Kohn-Sham 
equations of the usual form, because the resulting effective potential is dif- 
ferent for each orbital. As a consequence, various specialized algorithms for 
minimizing the PZ-SIC energy functional have been developed. For more 
details on these algorithms and some interesting applications in solid-state 
physics see Refs. [7TJ E2J [73] . For finite systems, PZ-SIC has also been im- 
plemented by means of the OEP [B3J [73], which produces a common local 
potential for all orbitals, and is discussed in the next paragraph. A detailed 
review of implementations and applications of PZ-SIC can be found in the 
contribution of Temmerman et al. in Ref. [TTJ. Alternatives to the PZ-SIC 
formulation of Ref. [93] have recently been analysed in [1181 H19j . with a 
view on either improving results obtained with PZ-SIC, or simplifying the 
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implementation of the correction. 

Since hybrid functionals, Meta-GGAs, SIC, the Fock term and all other 
orbital functionals depend on the density only implicitly, via the orbitals 
4>i[n], it is not possible to directly calculate the functional derivative v xc = 
5E xc /5n. Instead one must use indirect approaches to minimize E[n] and 
obtain v xc . In the case of the kinetic-energy functional T s [{</>j[n]}] this in- 
direct approach is simply the Kohn-Sham scheme, described in Sec. HJ In 
the case of orbital expressions for E xc the corresponding indirect scheme is 
known as the optimized effective potential (OEP) |120j or, equivalently, the 
optimized-potential model (OPM) [121j . The minimization of the orbital 
functional with respect to the density is achieved by repeated application of 
the chain rule for functional derivatives, 



54>i(r') Sv s (r") Sn(r 

(91) 

where E°^ b is the orbital functional (e.g., the Fock term) and v s the KS 
effective potential. Further evaluation of Eq. ( 19TT) gives rise to an integral 
equation that determines the w xc [n] belonging to the chosen orbital func- 
tional ^[{^[n]}] [T201 [T22] . As an alternative to solving the full OEP 
integral equation, Krieger, Li and Iafrate (KLI) have proposed a simple but 
surprisingly accurate approximation that greatly facilitates implementation 
of the OEP [LH]. 

The application of the OEP methodology to the Fock term (1561) . either 
with or without the KLI approximation, is also known as the exact-exchange 
method (EXX). The OEP-EXX equations have been solved for atoms |120| 
11211 1123] and solids [124L 1125] . with very encouraging results. Other orbital- 
dependent functionals that have been treated within the OEP scheme are 
the PZ self-interaction correction [6^1 [71] and the Colle-Salvetti functional 
[123J. A detailed review of the OEP and its KLI approximation is Ref. |122j . 

The high accuracy attained by complex orbital functionals implemented 
via the OEP, and the fact that it is easier to devise orbital functionals than 
explicit density functionals, makes the OEP concept attractive, but the com- 
putational cost of solving the OEP integral equation is a major drawback. 
However, this computational cost is significantly reduced by the KLI approx- 
imation [120] and other recently proposed simplifications [1261 11271 tl28j . In 
the context of the EXX method (i.e., using the Fock exchange term as orbital 
functional) the OEP is a viable way to proceed. For more complex orbital 
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functionals, additional simplifications may be necessary |120l I12(j| 11271 fl28j . 

A further reduction of computational complexity is achieved by not eval- 
uating the orbital functional self-consistently, via Eq. fl9Tl) . but only once, 
using the orbitals and densities of a converged self-consistent LDA or GGA 
calculation. This 'post-GGA' or 'post-LDA' strategy completely avoids the 
OEP and has been used both for hybrid functionals and Meta-GGAs [108, 
11091 11121 1113] . A drawback of post methods is that they provide only ap- 
proximations to the selfconsistent total energies, not to eigenvalues, effective 
potentials, orbitals or densities. 

In the case of hybrid functionals, still another mode of implementation has 
become popular. This alternative, which also avoids solution of Eq. (1911) . is 
to calculate the derivative of the hybrid functional with respect to the single- 
particle orbitals, and not with respect to the density as in (191|) . The resulting 
single-particle equation is of Hartree-Fock form, with a nonlocal potential, 
and with a weight factor in front of the Fock term. Strictly speaking, the 
orbital derivative is not what the HK theorem demands, but rather a Hartree- 
Fock like procedure, but in practice it is a convenient and successful approach. 
This scheme, in which self-consistency is obtained with respect to the single- 
particle orbitals, can be considered an evolution of the Hartree-Fock Kohn- 
Sham method [6], and is how hybrids are commonly implemented. Recently, 
it has also been used for Meta-GGAs [2]. For occupied orbitals, results 
obtained from orbital selfconsistency differ little from those obtained from 
the OEP. 

Apart from orbital functionals, which are implicit nonlocal density func- 
tionals because the orbitals depend on the density in a nonlocal way, there 
is also a class of explicit nonlocal density functionals. Such nonlocal density 
functionals take into account, at any point r, not only the density at that 
point, n(r), and its derivatives, Vn(r) etc., but also the behaviour of the 
density at different points r' ^ r, by means of integration over physically 
relevant regions of space. A typical example is 



where e!f c m is the per-particle xc energy of the homogeneous electron liquid 
(see footnote 47). In the LDA one would have n(r) = n(r), but in the 
average-density approximation (ADA) one takes [129J 




(92) 




(93) 
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where w [n] (| r — r ; | ) is a weight function that samples the density not only 
semilocally, as do the GGAs, but over a volume determined by the range of 
w. Conceptually similar to the ADA is the weighted-density approximation 
(WDA) [129] . In terms of the pair-correlation function (see Sees. 13.5.21 and 
14.1.21) the LDA, ADA and WDA functionals can be written as 



EL DA [n] = \\<?r\ dV (^[n(r)](r 
Ei DA [n] = jjd 3 rj dV (£w[n(r)](r 
EZ DA [n] = jjd*rj dV (g hom [n(v))(v - 



r') - 1) 


(94) 


r') - 1) 


(95) 


r') - 1) , 


(96) 



where in each case ghom( r — r ') is the pair-correlation function of the homo- 
geneous electron liquid, averaged over the coupling constant e 2 

The dependence of these functionals on n(r), the integral over n(r), in- 
stead of on derivatives, as in the GGAs, is the reason why such functionals 
are called nonlocal. In practice, this integral turns the functionals compu- 
tationally expensive, and in spite of their great promise they are much less 
used than GGAs. However, recent comparisons of ADA and WDA with LDA 
and GGAs for low- dimensional systems [1141 1130] and for bulk silicon |131] 
show that nonlocal integral-dependent density functionals can outperform 
local and semilocal approximations. 



6 Extensions of DFT: New frontiers and old 
problems 

Up to this point we have discussed DFT in terms of the charge (or par- 
ticle) density n(r) as fundamental variable. In order to reproduce the cor- 
rect charge density of the interacting system in the noninteracting (Kohn- 
Sham) system, one must apply to the latter the effective KS potential v s = 
v + vh + v xc , in which the last two terms simulate the effect of the electron- 
electron interaction on the charge density. This form of DFT, which is the 
one proposed originally [22], could also be called 'charge-only' DFT. It is not 
the most widely used DFT in practical applications. Much more common 
is a formulation that employs one density for each spin, nj(r) and n^(r), 
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i.e, works with two fundamental variables. In order to reproduce both of 
these in the noninteracting system one must now apply two effective poten- 
tials, v s ^(r) and f s ,].( r )0 This formulation of DFT is known as spin-DFT 
(SDFT) [89], ED]- Its fundamental variables rif(r) and nj.(r) can be used to 
calculate the charge density n(r) and the spin-magnetization density m(r) 
from 

n(r) = n-f(r) + nj.(r) (97) 
m(r) = /io(n T (r) -n|(r)), (98) 

where /io = qh/2mc is the Bohr magneton. More generally, the Hohenberg- 
Kohn theorem of SDFT states that in the presence of a magnetic field 
B(r) that couples only to the electron spin [via the familiar Zeeman term 
/ d 3 r m(r)B(r)] the ground-state wave function and all ground-state observ- 
ables are unique functionals of n and m or, equivalently, of nj and n^F^l 
Almost the entire further development of the HK theorem and the KS equa- 
tions can be immediately rephrased for SDFT, just by attaching a suitable 
spin index to the densities. For this reason we could afford the luxury of 
exclusively discussing 'charge-only' DFT in the preceding sections, without 
missing any essential aspects of SDFT. 

There are, however, some exceptions to this simple rule. One is the 
fourth statement of the HK theorem, as discussed in Sec. 13.21 Another is the 
construction of functionals. For the exchange energy it is known, e.g., that 
[152] 

E S x DFT [n h nj = \ [E° FT [2n,] + E^n,}) . (99) 

In analogy to the coordinate scaling of Eqs. ( J59l) - (!6Pj) . this property is of- 
ten called 'spin-scaling', and it can be used to construct an SDFT exchange 
functional from a given DFT exchange functional. In the context of the 
LSDA, von Barth and Hedin [89J wrote the exchange functional in terms 
of an interpolation between the unpolarized and fully polarized electron gas 
which by construction satisfies Eq. fl99l) . Alternative interpolation procedures 

52 More generally, one requires one effective potential for each density-like quantity to 
be reproduced in the KS system. Such potentials and corresponding densities are called 
conjugate variables. 

53 In the particular case B = the SDFT HK theorem still holds and continues to be 
useful, e.g., for systems with spontaneous polarization. In principle one could also use 
'charge-only' DFT to study such systems, but then n^(r) and n \ (r) become functionals of 
n(r) and nobody knows how to determine these functionals. 
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can be found in Ref. [H2]- GGA exchange functionals also satisfy Eq. 
by construction. For the correlation energy no scaling relation of the type 
(j99p holds, so that in practice correlation functionals are either directly con- 
structed in terms of the spin densities or written by using, without formal 
justification, the same interpolation already used in the exchange functional. 
In the case of the LSDA this latter procedure was introduced in Ref. 
and further analysed and improved in Ref. 
The Kohn-Sham equations of SDFT are 



(r) = e ia <j) ia (v), (100) 



where v sa (r) = v a (r) + ^(r) + v xc ^{r). In a nonrelativistic calculation the 
Hartree term does not depend on the spin label0 but in the presence of 
an externally applied magnetic field v a (r) = v(r) — o[i§B (where a = ±1). 
Finally, 



(101) 



In the presence of an internal magnetic field B xc (i.e., in spin-polarized sys- 
tems) v xc ^ — v xc ^ = fioB xc - This field is the origin of, e.g., ferromagnetism 
in transition metals. References to recent work with SDFT include almost 
all practical DFT calculation: SDFT is by far the most widely used form 
of DFT@ Some recent work on SDFT is described in Ref. |133j . A more 
detailed discussion of SDFT can be found in Refs. [SI El ED] , and a particu- 
larly clear exposition of the construction of xc functionals for SDFT is the 
contribution of Kurth and Perdew in Refs. [15l [16] . 

If the direction of the spins is not uniform in spaced one requires a for- 
mulation of SDFT in which the spin magnetization is not a scalar, as above, 
but a three-component vector m(r). Different proposals for extending SDFT 
to this situation are available |134l 1135} I136j . One mechanism that can give 
rise to noncollinear magnetism is spin-orbit coupling. This is another rela- 
tivistic effect [28], and as such it is not consistently treated in either DFT 



54 Spin-spin dipolar interactions are a relativistic effect of order (1/c) 2 , as are current- 
current interactions. 

55 SDFT has become synonymous with DFT to such an extent that often no distinction 
is made between the two, i.e., a calculation referred to as a DFT one is most of the time 
really an SDFT one. 

56 Such 'noncollinear magnetism' appears, e.g., as canted or helical spin configurations 
in rare-earth compounds, or as domain walls in ferromagnets. 
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or SDFT. A generalization of DFT that does account for spin-orbit coupling 
and other relativistic effects is relativistic DFT (RDFT) |137l 1138] . Here the 
fundamental variable is the relativistic four-component current j^. RDFT 
requires a more drastic reformulation of DFT than does SDFT. In particu- 
lar, the KS equation of RDFT is now of the form of the single-particle Dirac 
equation, instead of the Schrodinger equation. There are also many subtle 
questions involving renormalizability and the use of the variational princi- 
ple in the presence of negative energy states. For details on these problems 
and their eventual solution the reader is referred to the chapters by Engel 
et al. in Refs. [10] and [19], and to the book by Eschrig [18]. A didactical 
exposition of RDFT, together with representative applications in atomic and 
condensed- matter physics, can be found in the book by Strange [28], and a 
recent numerical implementation is presented in Ref. [139J. 

To study the magnetic properties of matter one would often like to be able 
to obtain information on the currents in the system and their coupling to pos- 
sible external magnetic fields. Important classes of experiments for which this 
information is relevant are nuclear magnetic resonance and the quantum Hall 
effects. SDFT does not provide explicit information on the currents. RDFT 
in principle does, but standard implementations of it are formulated in a 
spin-only version, which prohibits extraction of information on the currents. 
Furthermore, the formalism of RDFT is considerably more complicated than 
that of SDFT. In this situation the formulation of nonrelativistic current- 
DFT (CDFT), accomplished by Vignale and Rasolt [140] 1141] , was a major 
step forward. CDFT is formulated explicitly in terms of the (spin) density 
and the nonrelativistic paramagnetic current density vector j p (r). Some re- 
cent applications of CDFT are Refs. [HI EM CHS] . E. K. U. Gross and 
the author have shown that the existence of spin currents implies the exis- 
tence of a link between the xc functionals of SDFT and those of CDFT [146J. 
Conceptually, this link is similar to the one of Eq. (|99|) between functionals 
of DFT and SDFT, but the details are quite different. Some approximations 
for xc functionals of CDFT are discussed in Refs. [146[ I147[ 1148] . 

In addition to SDFT, RDFT and CDFT, there exist many other gener- 
alizations of DFT that were designed for one or other special purpose. As 
examples we mention superconductivity [149[ 11501 1151[ I152j and spin-density 
waves |136l I153j . but there are many more [5-19]. For reasons of space we 
cannot discuss these extensions here. Instead, let us take a brief look at a 
problem that requires more radical departures from the framework of con- 
ventional DFT: excited states. DFT is formulated in terms of ground-state 
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densities, and it is not immediately obvious how one could extract informa- 
tion on excited states from them (although at least in the case of 'charge-only' 
DFT the fourth substatement of the HK theorem guarantees that this must 
be possible). 

Apart from the ad hoc identification of the KS eigenvalues with true exci- 
tation energies, there exists a considerable variety of more sound approaches 
to excited states in DFT that have met with some success. The early sug- 
gestion of Gunnarsson and Lundqvist [90] to use a symmetry-dependent xc 
functional to calculate the lowest-energy excited state of each symmetry class 
has been implemented approximately by von Barth |154j . but suffers from 
lack of knowledge on the symmetry dependence of the functional. More 
recent work on this dependence is Ref. [155J. An alternative approach to 
excited states, not restricted to the lowest energy state of a given symmetry, 
is ensemble DFT, developed by Theophilou [26] and further elaborated by 
Oliveira, Gross and Kohn [156]. In this formalism the functional depends on 
the particular choice for the ensemble, and a simple approximation for this 
dependence is available [156] . Some applications of this method have been 
worked out by Nagy |157j . 

Other DFT approaches to excited states can be found in Refs. [158j . [159] . 
[160] and [31], but the most widely used method today is time-dependent 
DFT (TD-DFT). The time-dependent generalization of the HK theorem, the 
Runge-Gross theorem, cannot be proven along the lines of the original HK 
theorem, but requires a different approach [5TI I161j . For recent reviews of 
TD-DFT see Ref. P-62J. Excited states have first been extracted from TD- 
DFT in Refs. [163[ 1164] . This approach is now implemented in standard 
quantum-chemical DFT program packages [1651 1166] and is increasingly ap- 
plied also in solid-state physics [70]. Another important application of TD- 
DFT is to systems in external time-dependent fields, such as atoms in strong 
laser fields |167l 1168] . First steps towards studying dynamical magnetic phe- 
nomena with TD-SDFT have been taken very recently [169] . 

All these extensions of DFT to time-dependent, magnetic, relativistic and 
a multitude of other situations involve more complicated Hamiltonians than 
the basic ah initio many-electron Hamiltonian defined by Eqs. to <^j. 
Instead of attempting to achieve a more complete description of the many- 
body system under study by adding additional terms to the Hamiltonian, 
it is often advantageous to employ the opposite strategy, and reduce the 
complexity of the ah initio Hamiltonian by replacing it by simpler models, 
which focalize on specific aspects of the full many-body problem. Density- 
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functional theory can be applied to such model Hamiltonians, too, once a 
suitable density-like quantity has been identified as basic variable. Following 
pioneering work by Gunnarsson and Schdnhammer |170j . LDA-type approx- 
imations have, e.g., recently been formulated and exploited for the Hubbard 
[171] . the delta-interaction |172j and the Heisenberg [173] models. Common 
aspects and potential uses of DFT for model Hamiltonians are described in 

Still another way of using DFT, which does not depend directly on ap- 
proximate solution of Kohn-Sham equations, is the quantification and clari- 
fication of traditional chemical concepts, such as electronegativity [6], hard- 
ness, softness, Fukui functions, and other reactivity indices [SJ 1175] , or aro- 
maticity [176] . The true potential of DFT for this kind of investigation is 
only beginning to be explored, but holds much promise. 

All extensions of DFT face the same formal questions (e.g., simultaneous 
interacting and noninteracting f-representability of the densities, nonunique- 
ness of the KS potentials, meaning of the KS eigenvalues) and practical 
problems (e.g., how to efficiently solve the KS equations, how to construct 
accurate approximations to E xc , how to treat systems with very strong cor- 
relations) as do the more widely used formulations 'charge-only' DFT and 
SDFT. These questions and problems, however, have never stopped DFT 
from advancing, and at present DFT emerges as the method of choice for 
solving a wide variety of quantum mechanical problems in chemistry and 
physics — and in many situations, such as large and inhomogeneous sys- 
tems, it is the only applicable first-principles method at all. 

The future of DFT is bright [31 EH, 1177] - but to be able to contribute 
to it, the reader must now leave the present superficial overview behind, and 
turn to the more advanced treatments available in the literature [5-19]. 
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